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Foreword to Reports on Computational 
Physiology 


Dear Reader, 


In 2016, Springer and Simula launched an Open Access series called the Simula 
SpringerBriefs on Computing. This series aims to provide concise introductions 
to the research areas in which Simula specializes: scientific computing, software 
engineering, communication systems, machine learning and cybersecurity. These 
books are written for graduate students, researchers, professionals and others who 
are keenly interested in the science of computing. We know that entering a new field 
of research or getting up-to-date on a new topic of interest can be very demanding 
and time consuming, for students and experienced researchers alike. Each volume 
presents a compact, state-of-the-art disciplinary overview and raises essential critical 
questions in the field, all in approximately 100 pages. 


Simula’s focus on computational physiology has grown considerably over the last 
decade. Our researchers collaborate with partners around the world in interdisci- 
plinary teams to develop multi-scale mathematical models of excitable tissues (brain 
and heart). These models are becoming increasingly complex and accurate, in par- 
ticular as they are compared to experimental and clinical data. Since 2014, the 
University of California, San Diego (UCSD) and Simula have organized an annual 
summer school in computational physiology, in which graduate students spend the 
second two weeks of June in Oslo learning the principles underlying mathematical 
models commonly used to study the heart and the brain. The students are then as- 
signed a research project to work on over the summer. In August the students travel to 
the University of California, San Diego to present their findings. Each year, we have 
been impressed by the students’ abilities to learn the huge amount of mathematics 
and physiology theory required for their research projects, the results of which often 
contain the rudiments of a scientific paper. 


As a result of our expanding activities in this field, we have decided to publish a 
branch of the SimulaSpringer Briefs that is specifically focused on computational 
physiology. Each volume in this series will explore multiple physiological questions 
and the models developed to address them. Each of the questions will, in turn, be 
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packaged into a short report format (6-10 pages) that provides a succinct summary 
of the findings and, whenever possible, the software used will be made publicly 
available. All reports in this series are subjected to peer-review. We would like to 
emphasise that we do not require that reports represent new scientific results; rather, 
they can reproduce or supplement earlier computational studies or experimental 
findings. 


The main driver for this series is to enable the publication of project reports from 
the annual summer school in Computational Physiology, but we will also publish 
related reports that fit the overall purpose of the SimulaSpringer Briefs. Due to the 
Covid-19 pandemic, the summer school in 2020 had to be cancelled, and as such 
this first issue of the new series is not a collection of project reports. However, the 
topic of the first issue fits very well in the framework of computational physiology; 
it presents novel methods and software for the simulation of excitable cells. 


It is a pleasure to thank our collaborators at SpringerNature for their superbly efficient 
handling of this manuscript. In particular, we are grateful for the sound advice and 
support from Dr. Martin Peters, the Executive Editor for Mathematics, Computational 
Science and Engineering. We would also like to thank Dr. Henrik Nicolay Finsberg 
for his excellent technical support in this project. 


Fornebu, Norway Dr. Kimberly J McCabe 
September 2020 Dr. Rachel Thomas 
Dr. Andrew D McCulloch 

Dr. Aslak Tveito 


Preface 


Partial differential equations (PDEs) have proved to be immensely useful in mod- 
elling Nature; virtually all fields of science have their own equations, and every 
field of engineering is based on mathematical models formulated in terms of PDEs. 
This is astonishing given the fact that no model, but the very simplest ones, can 
be studied using analytical (paper and pencil) techniques. Numerical computations 
have proved tremendously useful in order to understand models formulated in terms 
of PDEs, and it can be argued that the computer was invented for the purpose of 
solving such equations. The computer is extremely well suited to perform the huge 
amounts of tedious and highly repetitive computations that earlier had to be com- 
pleted by humans. However, for a very long time, the computers typically available 
at research labs could solve only simple models. In the eighties, PDEs was almost 
always studied in 1 or 2 spatial dimensions; the 2D geometry was very simple and 
the model was most often linear and scalar. That level of computational complexity 
allowed analysis of qualitative properties of PDEs, but was insufficient for studying 
realistic models in Science and Engineering. 


Over the past 30 years, we have witnessed a tremendous development in computing 
power both in terms of hardware, solution algorithms and software. This development 
has paved the way for realism in modeling; geometrical structures can now be 
represented with high degree of accuracy and complex systems of PDEs is applied 
to model the complex dynamics under consideration. This has led to greatly improved 
understanding throughout many branches of Science and had led to the development 
of considerably more accurate tools in Engineering. 


Computational electrophysiology is a branch of Science that has benefited greatly 
from developments in computational analysis of PDEs. This development started 
70 years ago with the celebrated paper by Hodgkin and Huxley (see (6)) and was 
followed 10 years later by the first cardiac action potential model developed by 
Noble (13). Since these groundbreaking papers, there have been intense efforts 
to understand how excitable cells works based on modeling and computations. This 
development has moved in tandem with new experimental techniques providing more 
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accurate data necessary for parameterizing the increasingly complex models. A large 
number of membrane models (many hundreds) have been developed (see e.g. (11)), 
and these models have been used together with the monodomain or bidomain models 
(see e.g. (5][21)) to study the electrochemical waves underpinning the contraction of 
the cardiac muscle. Similarly, the Cable equation (see e.g. (15)) has been extensively 
used to understand propagation of electrochemical signals by neurons. 


Both the monodomain and the bidomain models of electrophysiology are derived 
based on homogenization of the cardiac tissue. In the resulting models, this means 
that the both the extracellular space, the cell membrane and the intracellular space 
are assumed to exist everywhere in the computational domain. Specifically, this 
means that the cardiac cell is not explicitly present in these models. This approach 
to modeling cardiac tissue enable analysis of phenomena on a relatively large length 
scale (mm), but is useless when it comes to study processes going on at a small 
length scale (um). In 1993, the bidomain model was solved by Trayanova et al. ((20)) 
using 257 computational nodes, and at that time this was the best that technology 
would allow, and using a homogenized, large scale, model made perfect sense. More 
recently, however, models based on about 30 million mesh points are used allowing 
a characteristic mesh length of about 50um which is about half the length of a 
human ventricular cell. Further refinement of the mesh used in the monodomain 
and bidomain models is not very useful since converged solutions are obtained at 
a quite coarse mesh (~0.3mm, see e.g. [3)). This means that technology now 
allows simulation at a shorter length scale than the classical models (monodomain 
and bidomain) are meant to represent; it is impossible to gain information at the um- 
level using these models. Specifically, the cell is impossible to represent explicitly 
in the classical models and that clearly limits their usefulness. 


Interesting phenomena in electrophysiology take place close to the cell membrane. 
But since the cell is not explicitly present in the classical models, it is very difficult, 
if at all possible, to use such models to get a grip on what is going on in the vicinity 
of excitable cells. And since mesh resolution already has reached the um-scale, it is 
clearly about time to introduce the cell as the building block in models of excitable 
tissue. In fact, this development started several years ago and has been pursued by 
many authors; see e.g. (7][10] [I4] [19] [16] [18] [171 DA] [T). Recently, we have followed 
up on these papers aiming at developing models, algorithms and software for cell- 
based representation of excitable tissue. We represent the extracellular (E) domain, 
the cell membrane (M) and the intracellular domain (I) explicitly, and therefore refer 
to it as the EMI-model. 


In the first paper (see (23)) using the EMI model, we compared the results of the EMI 
model with the results of the classical Cable equations. In particular, we assessed 
the magnitude of the error introduced in the Cable equation by ignoring the ephaptic 
effects (i.e. assuming that the extracellular potential is constant). Then the same 
model was used to assess the effect on the action potential of placing a microelec- 
trode array in the extracellular domain (see (2)). Furthermore, we have developed 
computational techniques for solving the EMI equations (see (22][8)) and used it to 
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study properties of the conduction velocity of electrochemical wave traversing the 
cardiac muscle during a heart-beat (see (9)). Quite recently, the EMI model has been 
extended to also account for ion concentrations in the entire computational domain 
using the electroneutral Kirchhoff-Nernst-Planck (KNP) model, and the resulting 
model is referred to as the KNP-EMI model (see (4)). This opens the possibility for 
analysing the effect of depression waves traversing cortical tissue. 


The main of advantage of the EMI-models is the possibility for the modeler to 
represent local properties of the cell and to study dynamics in the vicinity of the 
cells at the zm-scale. This opens vast possibilities for deeper understanding of the 
dynamics of collections of excitable cells. However, a main disadvantage is that 
the model is more complex and more computationally demanding that the common 
monodomain, bidomain, and cable equations. The purpose of this edition of the 
Simula SpringerBriefs on Computing is to provide succinct introductions to various 
aspects of the EMI-models, the solution algorithms and the software used to study 
these models. These models are not straightforward to implement and we therefore 
think it is useful to provide software for anyone interested in using the models. Note 
that it is specifically not our intention here to provide substantial new contributions 
to developments of models, algorithms or software, but rather to aid readers by 
providing easy and readable accounts of this material. 


In order to avoid misunderstanding, we would like to add that we do not think this 
is the end of the monodomain-, bidomain-, or the cable-model. These models have 
been extremely useful and in combinations with membrane models they basically 
represent our collective knowledge about how excitable cells work. Much work 
remains to be done with these equations and the models we suggest are far too 
computationally demanding to be a realistic alternative for full scale simulations of 
human organs. Also, again in order to avoid misunderstandings, homogenization is 
still with us in the EMI models; the EMI models takes the typical length scale form 
mm to um, but the atomic scale is still 10000 smaller, so the models we use in both 
E, M and I all represent averages. 


Fornebu, Norway Aslak Tveito 
September 2020 Kent-Andre Mardal 
Marie Rognes 
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Chapter 1 


Derivation of a Cell-Based Mathematical Model 
of Excitable Cells 


Karoline Horgmo Jeger! and Aslak Tveito'? 


Abstract Excitable cells are of vital importance in biology, and mathematical mod- 
els have contributed significantly to understand their basic mechanisms. However, 
classical models of excitable cells are based on severe assumptions that may limit 
the accuracy of the simulation results. Here, we derive a more detailed approach to 
modeling that has recently been applied to study the electrical properties of both 
neurons and cardiomyocytes. The model is derived from first principles and opens up 
possibilities for studying detailed properties of excitable cells. We refer to the model 
as the EMI model because both the extracellular space (E), the cell membrane (M) 
and the intracellular space (I) are explicitly represented in the model, in contrast to 
classical spatial models of excitable cells. Later chapters of the present text will focus 
on numerical methods and software for solving the model. Also, in the next chapter, 
the model will be extended to account for ionic concentrations in the intracellular 
and extracellular spaces. 


1.1 Introduction 


Mathematical modeling has a great potential for increasing our understanding of the 
physiological processes underlying the function of the body. For example, modeling 
of the electrical properties of excitable cells may provide insight into the complex 
electrical signaling involved in a number of important functions, like transfer of 
information through neurons and coordination of the pumping of the heart. A popular 
model of the conduction of electrical signals in neurons is the so-called cable equation 
(30), whereas the extracellular potential surrounding neurons is often modeled using 
the point-source or line-source approximations (6). The three aforementioned models 
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have been used extensively to gain insight into the function of neurons and the 
interpretation of measurements of the extracellular potential around neurons (6} 
[12). Correspondingly, the conduction of electrical signals through the heart is 
traditionally modeled using the homogenized bidomain and monodomain models 
(20). These models are also widely used and have, for instance, provided insight into 
mechanisms of cardiac arrhythmias (26] [33). 


However, despite the success of the above-mentioned classical models of excitable 
cells, the models have certain shortcomings that may make them inaccurate or im- 
practical in some situations. For example, in the derivation of the cable equation, the 
extracellular potential is often assumed to be constant (30][15). Therefore, changes in 
the extracellular potential generated by the neuron itself or by neighboring neurons 
(i.e., ephaptic effects) are ignored in the model. This could potentially lead to in- 
accuracies (16][2] [34). In addition, the point-source and line-source approximations 
rely on the assumption that the extracellular space is infinite and homogeneous. Con- 
sequently, the models might not be well-suited to interpret extracellular potentials 
measured when large measurement electrodes are present in the extracellular space 
close to the neurons (3). 


Moreover, the bidomain and monodomain models represent cardiac tissue in a ho- 
mogenized manner, assuming that the intracellular space, the extracellular space 
and the cell membrane exist everywhere in the tissue. Because the geometry of the 
individual cells is not represented, it is very hard to use the models to study the effect 
of, e.g., the cell geometry or a non-uniform distribution of ion channels on the cell 
membrane. These properties are both believed to influence cardiac conduction, but 
their exact effects are not fully understood and call for further investigations (28][21). 
In addition, it has been proposed that ephaptic coupling between cardiac cells might 
occur at small extracellular clefts located at the intercalated discs between cells (29). 
Since the geometry of the extracellular space is not represented in the homogenized 
models, it is difficult to use these models to study such ephaptic effects. 


In order to account for the difficulties related to the classical models, alternative 
electrophysiological models have been developed (e.g., (28][31][24)). In this chapter, 
we consider one of these alternative models, referred to as the EMI model, because 
it explicitly represents the extracellular space (E), the cell membrane (M) and the 
intracellular space (T). This model has been used to study both neurons 
and cardiac tissue (32][31][18). Because the model represents the extracellular space, 
the membrane and the intracellular space in a coupled manner, the model allows for 
representation of ephaptic effects between neurons or cardiomyocytes (18). In 
addition, since the geometry of the extracellular space is explicitly represented, the 
model allows for representation of non-homogeneous extracellular surroundings (3). 
Furthermore, since the geometry of each cell is represented, the model allows for 
study of the effect of cell geometry and non-uniform distributions of ion channels 
on cardiac conduction properties (18). 


In other words, the EMI model allows for a more detailed representation of excitable 
cells and tissues than classical models of computational electrophysiology. In fact, 
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Fig. 1.1: Illustration of an EMI model domain consisting of an extracellular domain, 
Qe, a cell membrane, I’, and an intracellular domain, €. 


the classical models mentioned above can be derived from the more detailed EMI 
model by introducing certain simplifying assumptions, see e.g., [1 1). In this 
chapter, however, we focus on deriving the EMI model from Maxwell’s equations of 
electromagnetism. 


1.2 Derivation of the EMI Model 


In this section, we present a derivation of the EMI model for excitable cells. This 
derivation is to a large extent based on the derivation found in (1}[17). We consider 
a domain separated into an extracellular part, Qe and an intracellular part, Q;, like 
illustrated in Figures 1.1 and 1.3. The cell membrane, denoted by I, is defined as the 
boundary between Q; and Qe. Here, we derive a model for the electrical potentials 
in both a domain with a single cell (Figure 1.1) and in a domain with two cells 
connected at an intercalated disc denoted by Tj » (Figure 1.3). 


1.2.1 Fundamental Equations 


We base the derivation of the EMI model on two of the quasi-static approximations 
of Maxwell's equations, i.e., 


VxE=0, (1.1) 
VxH-J. (1.2) 
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Here, E is the electric field (typically in uF/cm), H is the magnetic field (typically in 
uA/cm) and J is the density of free currents (typically in wA/cm7). In the quasi-static 
approximation of (1.2), it is assumed that free unbalanced charges are instantly bal- 
anced. The assumptions hold in the intracellular and extracellular spaces. However, 
we assume that charges may accumulate at the cell membrane and at the intercalated 
discs between cells. Therefore, we let (1.2) at these locations be replaced by the 
corresponding equation without the quasi-static approximation, i.e., by 


OE 
VxH- a 1. 
x Jtege (1.3) 


where & is the permittivity of the medium (typically in uF/cm). In addition, we 
assume that Ohm's law applies in the intracellular and extracellular domains. This 
means that 

J =cE, (1.4) 


where o is the conductivity of the considered medium (typically in mS/cm). We also 
note that (1.1) implies that E is a conservative vector field and that it therefore can 
be defined as the gradient of a scalar field (9). More specifically, we can define 


E = -Vu, (1.5) 


where the scalar u is the electric potential (typically in mV). 


1.2.2 Model for the Intracellular and Extracellular Domains 


In order to derive equations for the intracellular and extracellular domains, we take 
the divergence of both sides of (1.2) and apply the vector identity V - (V x H) = 0, 
which holds for any vector H (9). This yields 


V-J=0. 
Inserting (1.4) and (1.5), we obtain the Laplace equation 
V-oVu =0. (1.6) 
More specifically, for the intracellular and extracellular domains, we have 


V.o;Vuj-0 in Q,, (1.7) 
V-oeVue =0 in Qe, (1.8) 
where c; and ce are the intracellular and extracellular conductivities, respectively, 


and u; and ue are the electric potentials in the intracellular and extracellular domains, 
respectively. 
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1.2.3 Model for the Membrane 


In order to derive the EMI model equations for the membrane, we consider a volume 
element, B, intersected by the membrane. This volume element may be divided into 
an extracellular part, Be, and an intracellular part, B;, such that Be U B; = B and 
Be N B; = 0, as illustrated in Figure 1.2A. In each of these domains, we assume 
that (1.3) holds. Taking the divergence of both sides of (1.3) and applying the vector 
identity V - (V x H) = O results in 


V«ls-Vig—. (1.9) 


Integrating this equation over each of the volume elements B; and Be, we get 


f vJav-- f Meo ay 
Bi Bi ot 


l vJav--f Ue av, 
Be Be Ot 


and applying the divergence theorem (see e.g., (9)), we obtain 


OE 

J- ns, as=- f e— . ng, dS, (1.10) 
ðB; op, Ot 
OE 

Ings =- f e% . npg, dS. (1.11) 
ÓB. üB, Ot 


Here, ng, and ng, are the outward pointing normal vectors of B; and Be, respectively. 
Furthermore, in (1.10) and (1.11), the left-hand side terms represent the free ionic 
current and the right-hand side terms represent the capacitive current. 


1.2.3.1 Ionic Current 


We start by considering the left-hand side of (1.10), representing the ionic current. 
Here, we note that the boundary 0B; may be split into two parts, l'a and 0B; V Tg, 
where Ig is the part of ôB; coinciding with the membrane and 0B; \ Ip is the 
remaining part (see Figure 1.2A). We can then write 


| Jn ds = f Ja des f J- n; dS, (1.12) 
OB; OBi\TB Tg 


where nj is the outward pointing normal vector of the membrane and ng, is the out- 
ward pointing normal vector of the remaining part of B;, as illustrated in Figure 1.2A. 
At the membrane, the current density, J, consists of currents through different types 
of ion channels, pumps and exchangers located at the membrane. This current den- 
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Fig. 1.2: A: Illustration of a volume element, B, intersected by the membrane. The 
volume element is separated into an intracellular part, B;, and an extracellular part, 
Be. B: Illustration of a small volume element, B, located on the extracellular part of 
the membrane. 


sity is typically denoted by Tion and given in units of uA/cm?. By convention, Tion is 
defined to be positive for a flux of positive ions out of the cell (i.e., in the direction 
of nj). This gives 


Jonjds = f lion dS. (1.13) 
IN Tg 


The boundary 0B; V Tg is located in the intracellular domain. Here, we assume that 
the current density, J, is given by Ohm's law (1.4), such that 


f J- ng, as= [ c;E - ng, dS. (1.14) 
OB;NMp OB;NM'p 


Inserting (1.13) and (1.14) into (1.12), we get 
J: na, dS sf g;E - ng, dS +f Tion AS, (1.15) 
OB; OBi\'p Tp 


and similar arguments for the extracellular part of the membrane yield 


J-ng, dS i c, E - ng, dS -f lion dS. (1.16) 
OB, OB, Ng Tg 


Note that the negative sign in front of the last term is due to the fact that ne = —n;. 


1.2.3.2. Capacitive Current 


For the right-hand side part of (1.10), representing the capacitive current, we again 
split the integral into two parts 
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OE OE OE 
| €£—— ‘UB; as= f E> ‘UB; as+ f ep — n; ds. 
OB; Ot OBi\lB Ot Tp Ot 


Here, er is the permittivity of the membrane. Following the quasi-static assumptions, 


we assume that the term cf is negligible for the part of 0B; that does not coincide 
with the membrane. Furthermore, from (1.5), we get E- n; = -Vu nj « v/d, where 
V = uj — Ue (1.17) 


is the membrane potential and d is the thickness of the membrane. We assume that 
the membrane can be treated as a capacitor formed by two parallel plates separated by 
an insulator. In that case, the membrane capacitance per area is given by C,, = ep/d 
(13). Therefore, 


OE OE er OV Ov 
ES inp dtes muse eal 0535 08 
L'a "5s Ls a” a d at ir: a 99 


Similar arguments for the extracellular side yield 


OE Ov 
is dS =- m= d$, 1.1 
n ng, dS fe Ór S (1.19) 


where the change of sign again is due to the fact that ne = —nj. 


1.2.3.3 Collecting the Ionic and Capacitive Currents 


Collecting the ionic and capacitive currents by inserting (1.15)-(1.16) and (1.18)- 
(1.19) into (1.10)-(1.11), we obtain 


ð 
li cE np, as+ f hon dS =- f Cm Z d5, 
OB \Tp Tp Tp Ot 


a 
| cE- ng, dS -f Ton d$ = i Cas, 
OBe\TB Tg Tg ðt 


which can be rewritten to 


f cE ng, as=- f Im aS, (1.20) 
OB;N p Tg 


I TeE - ng, dS = | Ij, dS, (1.21) 
OBe\TB Tg 


where the total membrane current density Im is defined as 


ðv 


Im = Cm 
i: ôt 


Pho (1.22) 
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Fig. 1.3: Illustration of an EMI model domain consisting of two cells, Ql and a, 


connected at an intercalated disc, T1,2 and surrounded by an extracellular domain, 
Qe 


We now wish to rewrite (1.20)-(1.21) to a differential form. We note that we can 
divide any volume element, B, intersecting the membrane into a purely intracellular, a 
purely extracellular, and a membrane intersecting part. We also know that (1.7)-(1.8) 
hold in the purely intracellular and extracellular parts. Therefore, we are interested 
in equations (1.20)-(1.21) as the size of B approaches zero. For example, we may 
consider a small extracellular volume element shaped as a cylinder, as illustrated in 
Figure 1.2B. As the height, Ahg, of this cylinder approaches zero, the integral over 
OB, V Tg approaches the integral over l'a, and we therefore get 


| oE-np.ds~ | c, E - ng, dS. 
OBe\lB Tg 


Inserting this approximation into (1.21), we obtain 


1 c, E-ng,dS ad Im dS > TE: npg, = Im, 
Tg Tg 
and inserting (1.5) and ne = —ng,, we get 
Ge Vue : Ne = In. (1.23) 
Similar arguments for the intracellular part of the membrane yield 
— gi Vu; : hj = In, (1.24) 


where the negative sign is due to the negative sign in (1.20). Finally, combining 
(1.23) and (1.24), we obtain 


Ce Vue : ne = —O; Vu; - nj = In, (1.25) 


where Im = Cm 2v + Tion and v = uj — Ue (see (1.22) and (1.17)). 
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1.2.4 Model for the Intercalated Disc 


In some cases, we wish to model cells that are connected to each other, as illustrated 
in Figure 1.3. We then let the intercalated discs connecting the cells be represented 
as boundaries between the intracellular domains, like the membrane is a boundary 
between the intracellular and extracellular domains. Furthermore, we assume that the 
intercalated disc have capacitive properties like the membrane, and that gap junctions 
allow for currents between neighboring cells, in the same manner as ion channels 
allows for currents between the intracellular and extracellular spaces. Therefore, the 
derivation of equations for an intercalated disc follows the exact same lines as the 
derivation of the membrane equations. More precisely, for two connected cells, we 
define an intercalated disc potential, w, by 


w-ul-u (1.26) 


where u : and u? are the electric potentials in en and OF respectively. In addition, 
we define a total intercalated disc current density, 71,2, by 


2 = pud + Tap, (1.27) 
where Jgap is the current density through the gap junctions, with positive direction 
in the direction from ol to 0 C1,» is the capacitance of the intercalated disc, 
and C 152 is the capacitive current density of the intercalated disc. Furthermore, 
following the same arguments as for the derivation of the membrane equations, we 
end up with an analogue to (1.25) of the form 


o? Vu? -n = -ojVui -nl = Do, (1.28) 


representing the total current density across the interface. 


1.2.5 Models of the Ionic Currents 


Mathematical models of the ionic currents governing the membrane potential of 
excitable cells come in a large variety of versions; see for several hundred 
examples. The simplest possible model is just a passive current of the form Jion = 
const - v, followed by a third order polynomial model. More realistic models tend to 
be more complex and are usually written on the form 


N 
los = V i (1.29) 


i=l 
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given in uA/cm?. Here, the individual currents can usually be written on the form 
I; = I;(v,s), where v denotes the membrane potential, given by u;i — ue, and s 
denotes gating variables and ionic concentrations. The celebrated model of the 
action potential of a neuron presented by Hodgkin and Huxley (see (14)) can be 
written on this form, and so can the first model of a cardiac cell presented by Nobel 
(25). A comprehensive and readable introduction to models of the membrane ionic 
currents is given in the survey (27). 


Correspondingly, the ionic currents through gap junctions between neighboring 
cells are often modeled by a simple passive model of the form Iga, = const - w. 
More detailed models of voltage-dependent gap junction dynamics have also been 


introduced (see e.g., (101[35)). 


1.2.6 Summary of the Model Equations 


In summary, the EMI model for a single cell surrounded by an extracellular domain 
(as illustrated in Figure 1.1) is given by the equations (1.7), (1.8), (1.17), (1.22) and 
(1.25), that is 


V-a;Vu; 20 in Qj, (1.30) 
V-aeVue =0 in Qe, (1.31) 
CeVue:ne ——c;Vuj:nj = Im atl, (1.32) 
V =U; — Ue at T, (1.33) 

Óv 1 
ET = C, Um — lion) at I, (1.34) 


where u;, Ue and v are the intracellular, extracellular and membrane potentials, 
respectively, typically given in mV. Moreover, c; and ce are the intracellular and 
extracellular conductivities, respectively (typically in mS/cm), Cm is the membrane 
capacitance (typically in u/F/cm?), and T denotes the cell membrane. The ionic 
currents through channels, pumps and exchangers at the membrane are denoted by 
Lion and typically given in uA/cm?. 


If several cells are connected at intercalated discs, as illustrated for two cells in 


Figure [1.3] the system of equations must be extended to include equations for the 
currents between cells. For two cells, this extension consists of the equations 


c;Vul-n?-2-o;Vul:n]s Ii» ati (1.35) 


ul-ul-w at 1,5, (1.36) 


l 


1 
w= c da = lga) — atl, (1.37) 
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where, as above, T'1? is the intercalated disc, n! is the outward pointing normal 
vector of ol, n? is the outward pointing normal vector of oF and ul and u? are the 
intracellular potentials (typically in mV) of Q! and Q?, respectively. Furthermore, 
C,» is the specific capacitance of the intercalated disc (typically in uF/cm?), and 
Tgap is the current through the gap junctions (typically in pA/cm?). 


1.3 Conclusion 


In the present chapter, we have derived the EMI model. The EMI model predicts 
electrical potentials in cells with an explicit geometrical representation and thus 
allows for more detail than homogenized models of excitable tissue. In the next 
chapter (7, Chapter 2), the model will be extended by taking ion concentration in the 
extracellular and intracellular spaces into account. Numerical solutions of the EMI 
models will be presented in Chapter 4), Chapter 5) and Chapter 6). In 
these chapters the readers will also be pointed to open software that can be used to 
solve the EMI model. 


Open Access This chapter is licensed under the terms of the Creative Com- 
mons Attribution 4.0 International License (http://creativecommons.org/licenses/ 
by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in 
any medium or format, as long as you give appropriate credit to the original aut- 
hor(s) and the source, provide a link to the Creative Commons license and indicate 
if changes were made. 

The images or other third party material in this chapter are included in the chap- 
ter's Creative Commons license, unless indicated otherwise in a credit line to the 
material. If material is not included in the chapter's Creative Commons license and 
your intended use is not permitted by statutory regulation or exceeds the permitted 
use, you will need to obtain permission directly from the copyright holder. 
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Chapter 2 


A Cell-Based Model for Ionic Electrodiffusion in 
Excitable Tissue 


Ada J. Ellingsrud!, Cécile Daversin-Catty! and Marie E. Rognes! 


Abstract This chapter presents the KNP-EMI model describing ion concentrations 
and electrodiffusion in excitable tissue. The KNP-EMI model extends on the EMI 
model by removing the assumption that ion concentrations are constant in time 
and space, and may as such be more appropriate in connection with modelling 
e.g. spreading depression, stroke and epilepsy. The KNP-EMI model defines a system 
of time-dependent, nonlinear, mixed dimensional partial differential equations. We 
here detail the derivation of the system and present a numerical example illustrating 
how ion concentrations evolve during neuronal activity. 


2.1 Introduction and Motivation 


In this chapter, we present an extension of the EMI model, presented in Chapter 
1), describing ion concentrations and electrodiffusion in excitable tissue. The EMI 
model is based on the assumption that intra- and extracellular ion concentrations are 
constant in time and space. This is often a good approximation, as ion concentrations 
in healthy tissue typically quickly return to base levels after neuronal activity due to 
cellular mechanisms such as e.g. membrane pumps and glial cell buffering. However, 
there are scenarios where this assumption is inadequate. 


Several cerebral pathologies are associated with increased neuronal activity (3), such 
as e.g. seizures and epilepsy (10} 6} [1), stroke (17), and spreading depression (22). 
In particular, periods of neuronal hyperactivity can lead to substantial variations in 
extracellular ion concentrations. These variations will in turn (i) influence mem- 
brane reversal potentials and (ii) generate diffusive currents. Changes in the reversal 
potentials, caused by local ionic shifts, may affect the dynamical properties of the 
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neurons [24). On the other hand, diffusive currents, driven by ion concen- 
tration gradients, can shift the extracellular potential (8] B). Mathematical models 
addressing the aforementioned phenomena and pathologies should therefore also 
account for ion concentrations, their spatial and temporal gradients and associated 
dynamics. 


In this chapter, we derive a system of time-dependent, nonlinear partial differen- 
tial equations describing the distribution and evolution of ion concentrations in a 
geometrically-explicit representation of the intra- and extracellular domains using 
the electroneutral Kirchhoff-Nernst-Planck (KNP) model (21). We will refer to this 
model as the KNP-EMI model, see also e.g. (5). 


2.2 Derivation of the Equations 


Let the computational domain Q and subdomains €, Qe, and I be defined as in 
the previous chapter 1.1. For simplicity and clarity, we present the mathematical 
model for one intracellular region Q; = Q; with membrane I below. We model 
a set K of intracellular and extracellular ion concentrations, and note that key ions 
in excitable tissue are potassium (K^), sodium (Na*), and chloride (CI). For each 
ion species k € K and each region r € {i,e}, we model the ion concentrations 
e : Q, x (0,T] — R (mol/m?), and electrical potentials uy : Q, x (0,T] —> R 
(V), and additionally the total transmembrane current density Im : T x (0,T] —D R 
(A/m?). 


2.2.1 Equations in the Intracellular and Extracellular Volumes 


In the EMI model, the free current densities J;, Je (uA/cm?), c.f. (1.4), are assumed 
to satisfy Ohm's law. To include diffusive ion effects, we instead assume that the 
free current density is composed of flux density contributions J k (mol/(m?s)) from 
different ions K as: 

= >) FoR ino, (2.1) 


keK 


where z^ is the valence of ion species k and F (C/mol) is Faraday's constant. 
Furthermore, we assume that ions can move by diffusion and/or in response to the 
electrical field as charged particles. Hence, the ion flux densities are modelled as 
the sum of two terms: (i) the ion concentrations that are transported via electrical 
potential gradients o*Vu, and (ii) the diffusive movement of ions due to ionic 
gradients Dk V ck: 

JE =-o Vu, —D*VcF in Q,, (2.2) 


y 
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where D* (m?/s) and c denote the effective diffusion coefficient and the conduc- 

tivity for ion species k in region r, respectively. The conductivity o* depends on 

the concentration of ion species k and the diffusion coefficient D^ in the following 

manner: "un 
A Dez 

e; Sag) = 4 


in Q,. (2.3) 


r 


Here, the constant y = RT F^! combines Faraday's constant F, the absolute temper- 
ature T (K), and the gas constant R (J/(K mol)). Moreover, the bulk conductivity c. 
can be expressed as: 


F 
c, = o, (c*) = — >» DEck(z ino, (2.4) 
d keK 


See e.g. for a derivation of the conductivity (2.3) and the bulk conductivity (2.4). 
Comparing with (1.4) and (1.5), we note the dependency on the ion concentrations 
in the conductivity o, in (2.3), and the second term accounting for ion diffusion 
in (2.2). 


As in Chapter 1, we stipulate that: 


V.J; 20 in Qi, (2.5) 
V.Je 20 inQe. (2.6) 


Finally, conservation of ions for the bulk of each region Q, gives that: 


O[k]i 


Adi y. jk = in Q; 
3; tV =0 in, (2.7) 
a +V-Jk=0  inQ, (2.8) 


for t € (0,T]. 


2.2.2 Membrane Currents 


We next turn to modelling the cell membrane currents and membrane potential 
across the interface I'. As in Chapter 1, we introduce the membrane potential v as 
the jump in the electrical potential over the membrane: 


y = Ui — Ue onT. (2.9) 


We also introduce the total membrane current as the combination of a capacitive 
current and ion specific currents: 
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ov 
In = cap + Toon = Cm a, + lion. (2.10) 


where the total channel current Jion is the sum of the ion specific channel currents 
I: 


ion’ 


lac da Tp = we a (2.11) 
keK 


The channel currents rm are subject to modelling, and will be discussed briefly in 
Section 2.2.2.1. 


Using our concepts, we have that the total ionic current density Im : T x (0, T) — R 
(A/m?) across the interface I (from the intracellular to the extracellular domain) is 
given by: 
-F Xon =F YD) AS ny = Im. (2.12) 
keK keK 
It now remains to specify a set of interface conditions for the specific ion fluxes 
Jk - n, forr € (i, e). 


Here, we propose a heuristic approach via ion specific capacitive current modelling, 
and note that an alternative approach is presented in (15). As for the total current, 
we assume that the capacitive current can be represented as a sum of ion specific 


contributions: 
DNS (2.13) 
keK 


Without loss of generality, we let the ion specific capacitive current is , in region 


Q, at the interface [ be some fraction ak of the total capacitive current cap: 
I axe ak leap: (2.14) 
Specifically, we assume that: 


uic Dk(z*)?[k], . 
Diex DLHY 


and note that X zeg ak = ] forr € (i, e). By the above definitions, (2.10) and (2.12), 
we let the intracellular and extracellular ion fluxes across the membrane be given by: 


(2.15) 


y — P. T af (Im = lion) E ds I + ak (Im = lion) 
i i ; » ? PE 


Fk 3 (2.16) 


for k € K. 
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2.2.2.1 Modelling Specific Ion Channels 


The membrane channel currents 1k (v) for each ion species k are subject to mod- 


elling. These currents are typically expressed on the form: 
Fon) = giU - E“), (2.17) 


where pt is the conductivity, and E* is the ion specific reversal potential (or Nernst 
potential), given by: 


sk, (2.18) 


This Nernst potential depends on the concentration ratio, whereas the Nernst potential 
in models without explicit modelling of ion concentrations is constant. Typical 
models include synaptic input currents, passive neuronal leak channels, or e.g. the 
Hodgkin-Huxley model (9). For more details on membrane current models and 
modelling, see e.g. (18). 


2.2.3 Summary of KNP-EMI Equations 


The KNP-EMI model equations follow from inserting (2.1) into (2.5)-(2.6), com- 
bined with (2.7), (2.8), (2.9), (2.10), and (2.16), and read as follows. 


For each ion species k € K and each region r € {i,e}, find the ion concentrations 
ck : Q, x (0,T] — R (mol/m?), the electrical potentials u, : Q, x (0,T] > R (V), 


r 


and the total transmembrane current density Im : T x (0,T] — R (A/m?) such that!: 


V-(F 5 D =0 in Q,, (2.19) 
k 
ück 
2 +V-JE=0 in Q,, (2.20) 
-F i WE -ne = DIa ny = Im at T, (2.21) 
k k 
V = Ui — Ue at I, (2.22) 
Ov 1 
gr = gum o) at T, (2.23) 


where the ion flux density J‘ is given by (2.2), and Jion is subject to modelling. A 
set of initial and boundary or compatibility conditions will close the system. 


! Note that the additional negative signs in (2.19) and (2.21), compared with the corresponding 
equations in Chapter 1, result from our physically consistent definition of the ion flux density J* as 
the negative gradient, cf. (2.2). 
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2.3 Numerical Solution of the KNP-EMI Equations 


The KNP-EMI model defines a complicated system of time-dependent, nonlinear, 
mixed dimensional partial differential equations. The number of unknowns depends 
on the number of ion species modelled. Some of the variables exist in the intracellular 
and extracellular domains, while others live on the lower-dimensional membrane. 
This setting is numerically challenging and calls for advanced techniques. 


To solve the KNP-EMI model numerically, one may consider a finite difference 
scheme to approximate the time derivatives, a linearization of ion flux densities J : 
and fractions o£, a splitting scheme to handle active ion channel current models, and 
a finite element discretization in space. Such a solution algorithm is detailed in (5), 
and we refer the reader to this description for further details. 


2.4 Comparing KNP-EMI and EMI during Neuronal 
Hyperactivity 


Neurons are negatively charged relative to their environment, with a resting mem- 
brane potential of about —70 mV. This resting potential is maintained by low con- 
centrations of sodium ions (Na*) and high levels of potassium ions (K*) inside 
the cell (23). Action potentials (neuronal activity) are generated by the opening of 
sodium and potassium channels in the cell membranes. The ionic gradient will drive 
sodium into the cell and depolarize the cell membrane. Next, the potassium channels 
open causing an outflux of potassium which in turn repolarizes the cell. 


As a result, there is a continuous need to pump potassium into the intracellular 
space and sodium out to the extracellular space to restore the electrochemical gra- 
dient across the cell membrane. One of the key mechanisms for this process is the 
Na/K/ATPase pump. The Na/K/ATPase pump actively transports 3 Na* ions out of 
the cell and 2 K* ions into the cell [20). Several pathologies are associated 
with increased neuronal activity, e.g. seizures and epilepsy (10][6] [1), and spreading 
depression (22). In periods of neuronal hyperactivity, the Na/K/ATPase pumps may 
not be able to restore the concentrations to baseline levels. Consequently, the elec- 
trochemical gradients may be reduced, and silenced neuronal activity and cellular 
swelling may occur (13). 


The ion concentration gradients observed during neuronal hyperactivity thus yields 
a suitable setting for illustrating differences between the KNP-EMI and the EMI 
frameworks. In particular, we compare the two frameworks both during normal 
neuronal activity (firing rate of 1 Hz) and during hyperactivity (firing rate of 50 Hz). 
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2.4.1 Model Parameters and Membrane Mechanisms 


We consider two idealized axons, represented by two parallel, rectangular domains, 
surrounded by extracellular space in three dimensions. The diameter of each axon is 
2.0- 1077 m, and they are separated by 1.0- 1077 m of extracellular space. Parameter 
values are as listed in Table 2.1. We refer to the supplementary code for a complete 
description of the model set-up (4). 


KNP-EMI membrane mechanisms The membrane mechanisms in the KNP-EMI 
model, cf. (2.11), are modelled using the standard Hodgkin-Huxley model (9) com- 
bined with a model for the Na/K/ATPase pump (12), the KCC2 cotransporter (24) 
and the NKCCI cotransporter (24). The Na/K/ATPase pump current Tarp (A/m?) is 
modelled as: 


^ 


I 
(+ E + Be 


Inrp = Isrel, ce) = (2.24) 


where Î is the maximum pump strength and mg and my, denote the pump thresh- 
old for extracellular potassium and intracellular sodium, respectively. Further, the 
transmembrane currents generated by the KCC2 cotransporter /kcco (A/m?) and the 
NKCC! cotransporter /wkcci (A/m?) are modelled as: 


EU 
Ikcc2 = Skcc2 In^ K aR (2.25) 
Ce Ce 
1 KCI NacCl 
IxKcci = SNKCCI gag n n( Ko) IM ea x xS ». (2.26) 


e 


where Skcc? and Snxcc are the maximal cotransporter strengths. Moreover, the cell 
is stimulated by prescribing a synaptic input 7,,, of the form: 


1-09 
IE. = 8synHe = (v - EF), Q.27) 
where a (s) is the synaptic time constant, H is the Heaviside function for a small 
region on the left side of the axons, and gsyn = 1.25 - 107? S/m?. In summary, the 
membrane channel currents for sodium, potassium and chloride are modelled as: 


Tg. ck) = = ent (v = ENS + gm h(v = ENS + 37 arp + INKCCI + i 


IE (v, ck) = g (v — ES) + gEn*(v — ES) - 2Inrp + Inxcci + Ikcco 
FEA, Cr) = giu (v — EC) - 2Ixcci - Ikcco, 

where, gk ak and 2* is the leak conductivity and the maximal conductivity for ion 
species k, respectively, the Nernst potential E^ for ion species k is as described in 
Section 2.2.2.1, and the gating variables m,h and n are described by the standard 
Hodgkin-Huxley ODEs, see e.g. for details. 
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EMI membrane mechanisms For the EMI model, we apply the standard Hodgkin- 
Huxley model and stimulate the cell by prescribing an input current of the form (2.27); 
thus, the membrane channels currents are modelled as: 


ion) = giak — EN*) + gia (V — EX) + gial — EC) 


+ gm? h(v — EN) + gFn*(v — EF) + [arp + Lyn, 


where EX, EN? and EC! are calculated by (2.18) with the initial values from the KNP- 
EMI model for the sodium and potassium concentrations. Similarly, the bulk conduc- 
tivities c; and c are calculated by (2.4), and the net current from the Na/K/ATPase 
pump /ATp is given by (2.24). Finally, there is no contribution from KCC2 and 
NKCC1, as both cotransporters mediate ion transport without any net charge move- 
ment across the membrane. 


Parameter Symbol Value Unit Reference 
gas constant R 8.314 J/(K mol) 
temperature T 300 K (23| 
Faraday’s constant F 9.648 - 10^ C/mol (23) 
membrane capacitance Cin 0.02 F/m (24| 
Na* diffusion coefficient DN 1.33 - 107? m?/s (23) 
K* diffusion coefficient DE 1.96 - 107? m?/s 
CI diffusion coefficient pa 2.03 - 107? m?/s (23| 
intracellular immobile anions c? 110 mM 
extracellular immobile anions c? 10 mM 

valence of immobile anions ZA -1 

Na* leak conductivity gy 0.281 S/m? * 
K* leak conductivity gr 0.43 S/m? F 
CI- leak conductivity gt 0.2 S/m? * 
K* HH max conductivity ak 360 S/m? (9) 
Nat HH max conductivity — g"* 1200 S/n? (9) 
maximum pump strength l 0.18 A/m? 
maximum KCC2 strength Sxcc2 0.0034 A/m? * 
maximum NKCCI strength Synxcci 0.023 Alm? $ 
ECS K* pump threshold my 3 mM * 
ICS Na* pump threshold Mat 12 mM * 
synaptic time constant a 1.0-107 s 

global time step At 10-1075 s 

local time step At* At/25 s 

spatial resolution Ax = Ay 2.5.1077 m 


Table 2.1: The physical and model parameters used in the simulations. The values 
are collected from Sterratt et al. (23), Hodgkin et al. (9), Wei et al. (24), whereas the 
values marked with * are computed by a steady state estimation. 
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The initial conditions for the intra- and extracellular ion concentrations, the mem- 
brane potential and the gating variables are listed in Table At the exterior 
boundary, we apply no flux boundary conditions for each ion species. 


Parameter Symbol Value Unit Reference 
initial intracellular Na* concentration EET: mM 

initial extracellular Na* concentration ceo 120 mM 

initial intracellular K* concentration cE? 80 mM 

initial extracellular K* concentration cK? 4 mM 

initial intracellular Cl* concentration p 7 mM 

initial extracellular Cl* concentration coh 112 mM 

initial membrane potential y? -67.74 - 10? V 

0 ay v?) 


initial HH gating value (Na* activation) m a R5 7 


anv^) E 
aq (v9 By (v?) 
0 av^) = 
as (v9) B, (v9) 


initial HH gating value (Na* inactivation h? 


Ea mE 


initial HH gating value (K* activation) n 


Table 2.2: Initial conditions. The initial ion concentrations are chosen such that the 
Nernst potentials are equal to those in the Hodgkin-Huxley model (9). The membrane 
potential is computed by a steady state estimation. 


2.4.2 Results and Discussion 


During normal activity, the KNP-EMI and the EMI models behave similarly, both 
for the membrane potential and the extracellular potential (Figure 2.1 A, B). The 
stimuli current depolarizes the membrane potential above the threshold for firing, 
and an action potential is initiated (Figure 2.1 A). Simultaneously, the extracellular 
potential decreases by ~ 0.13 mV, before quickly returning to baseline (Figure 2.1 
B). 


During hyperactivity, the KNP-EMI and EMI models differ (Figure 2.1 C, D, E, 
F). In both models, repeated action potentials are triggered. But, for the KNP-EMI 
model, we observe changes in the membrane potential between hyperpolarization 
phases. In particular, we conclude that the KNP-EMI membrane resting potential 
increases with repeated firing: after 5 action potentials (at t = 90 ms) the membrane 
potential has a minimum value of —75 mV, which is an 9% increase from the first 
action potential. Eventually, the membrane is depolarized to the point where action 
potentials can long longer be fired (Figure 2.1 E). 


The observed changes are caused by alterations in the ion concentration gradients. 
For each action potential, the extracellular Na* concentration decreases by 0.15 mM 
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and the extracellular K* concentration increases by 0.16 mM (Figure 2.2 A, B). Dur- 
ing normal activity (Figure 2.2 A, B), the ion concentrations will slowly be pumped 
back toward baseline levels, and the membrane potentials are not substantially af- 
fected by the small ion concentration changes. However, in the case of hyperactivity, 
the membrane mechanisms (i.e. pumps and cotransporters) are not able to keep up. 
Consequently, the extracellular Na* concentration will keep decreasing and the ex- 
tracellular potassium will keep increasing, causing the cell to depolarize (Figure 2.2 
C, D). 


In the KNP-EMI model (Figure 2.2 A, B), we note that 7.92 % of the extracellular 
K* concentration is restored, and 7.3 % of the extracellular Nat concentration is 
restored after 100 ms. That is, the extracellular concentrations do not reach baseline 
levels within the simulation period. Other studies have reported that it takes on the 
order of minutes (0.5 minutes (19), 6 minutes (2)) before the concentrations return 
to baseline after neuronal activity. 


2.5 Conclusions and Outlook 


In this chapter, we have presented a mathematical model, the KNP-EMI model, 
for ionic electrodiffusion in excitable tissue with an explicit representation of the 
intracellular, extracellular and membrane domains. For further reading on method- 
ological aspects, we refer to and references therein. This model extends on the 
EMI model presented in Chapter 1 and may be more accurate in situations with rapid 
and persistent changes in ion concentrations. Moreover, the KNP-EMI framework 
allows for modelling ligand-gated ion channels (e.g. NMDA receptors). 


The complexity of the KNP-EMI system yields a number of numerical challenges. 
The mere number of unknowns result in large systems of equations calling for 
efficient solution techniques. The nonlinearities in the system can easily lead to 
non-convergence and thus call for robust algorithms. Moreover, the coupling of full 
and lower dimensional domains and fields calls for well-posed numerical methods 
together with suitable simulation software. Further, the system couples different time 
scales: from neuronal action potentials taking place at the microscale to the slower 
diffusion process. In short, modelling ionic electrodiffusion in the EMI setting is an 
area with vast opportunities for further research. 
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Fig. 2.1: Comparison of potentials over time at fixed points in space predicted by 
the KNP-EMI and the EMI frameworks during normal activity (upper panels) and 
during hyperactivity (mid and lower panels). The membrane potentials for KNP-EMI 
and EMI during normal activity (A) and hyperactivity (C, E, F), and the extracellular 
potentials for KNP-EMI and EMI during normal activity (B) and hyperactivity (D). 
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Fig. 2.2: Time development of extracellular ion concentrations at a fixed point 
in space for the KNP-EMI framework during normal activity (upper panels) and 
hyperactivity (lower panels). The extracellular sodium (A) and potassium (B) con- 
centrations during normal activity, and the extracellular sodium (C) and potassium 
(D) concentrations during hyperactivity. 
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Chapter 3 


Modeling Cardiac Mechanics on a Sub-Cellular 
Scale 


Ashild Telle!, Samuel T. Wall! and Joakim Sundnes! 


Abstract We aim to extend existing models of single-cell mechanics to the EMI 
framework, to define spatially resolved mechanical models of cardiac myocytes 
embedded in a passive extracellular space. The models introduced here will be 
pure mechanics models employing fairly simple constitutive laws for active and 
passive mechanics. Future extensions of the models may include a coupling to the 
electrophysiology and electro-diffusion models described in the other chapters, to 
study the impact of spatially heterogeneous ion concentrations on the cell and tissue 
mechanics. 


3.1 Introduction 


A vast range of models have been developed for the force development of cardiac 
and skeletal muscle, on the scale of a single cross bridge (10), myofilament (3), sar- 
comere (2), and the complete cell (13). The scales involved and the main functional 
units considered on each scale are schematically illustrated in Figure 3.1. Common 
to most existing models is the fact that they focus on a single spatial scale, and any 
coupling between scales is fairly crudely represented. As an example, the model by 
Rice et al. is essentially a model of a single sarcomere (Fig. 3.1 D), which is 
normalized and then scaled to yield a realistic force output for cell- and tissue-level 
mechanics applications. Other models provide detailed descriptions of mechanisms 
and interactions on a molecular level (Fig. 3.1 F)(4][3), and are able to capture many 
of the characteristic non-linearities of muscle cell mechanics. However, key aspects 
of mechanical activation and force-length relationships are still not fully understood, 
and they may be the result of interactions between individual sarcomeres and myofib- 
ril bundles. A few attempts have been made at modeling interactions at this scale, and 
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A 


10 nm 


Fig. 3.1: The heart (A) is mainly composed of cardiac muscle cells, also called 
cardiomyocytes (B). Each cell (C) is composed of long tubes of sarcomeres (D), in 
which the thin and thick myofilaments overlap in layers (E). The interaction between 
these (F) causes the cardiac muscle to contract in a process called the cross-bridge 
cycle. 


have shown potentially interesting emergent behaviours (21 [11). Furthermore, heart 
failure and other pathologies are linked with heterogeneous intracellular calcium 
concentration resulting from disruptions of the calcium regulation system. Describ- 
ing the effect of such heterogeneities on the cell contraction and force development 
requires spatially resolved mechanics models on the sub-cellular scale. 


Finite element models of contracting myocytes have been proposed (8}[14), and have 
been used to explore the impact of model assumptions, calcium heterogeneity, and 
boundary conditions. The model presented by Ruiz-Baier et al. describes the 
individual myocyte as a hyperelastic material, and uses an active strain approach 
to describe the contraction. Both the passive and active mechanical properties are 
assumed to be homogeneous, but sub-cellular heterogeneities can easily be intro- 
duced. We here propose to extend the single myocyte model in to include the 
extracellular domain, and to model collections of cells, based on similar ideas used 
for the electrophysiology model presented in and (7] Chapter 1). 
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3.2 Models and Methods 


The motion and deformation of the heart can be described by the classical theory of 
non-linear solid mechanics. The primary unknown in our computational model will 
be the displacement vector u, which for each material point describes the difference 
between its current and original position. We have u = x — X, where X is the original 
(reference) position of a point, and x is its position after the deformation. From 
the displacement vector we can define the deformation gradient F = 0x/0X = 
I + Qu/OX, which is an essential quantity describing the deformation of a solid. See 
for instance (6) for a detailed introduction to non-linear solid mechanics. 


A characteristic feature of the heart and other muscles is that they contract and deform 
even in the absence of external loads. The overall deformation and mechanical state 
of the heart depends both on this active contraction and on the passive mechanical 
properties of the tissue. There are two main approaches for modeling the coupling 
of active and passive mechanics in cardiac tissue, often referred to as active strain 
and active stress. Both approaches are based on modeling the active and passive 
contributions separately, then combining them into a complete coupled model. 


In the active strain approach, the active-passive coupling is incorporated through a 
multiplicative decomposition of the deformation gradient F into active and passive 
components, F = Fp Fa. Here, Fa represents an active deformation governed by the 
cell state, and F, is a passive elastic deformation which ensures compatibility with 
loads and kinematic boundary conditions. The active stress approach is based on an 
additive split of the stress tensor into its active and passive components. In terms of 
the first Piola-Kirchhoff stress tensor P, the stress is written as P = Pp + Pa, where 
P, is a function of the cellular activation state and P, is a standard elastic stress 
derived from a strain energy function. 


Both of these approaches have their strengths and weaknesses. In general, the active 
strain approach is considered to be more suitable for deriving mathematically well- 
behaved constitutive laws, while the active stress concept is more easily coupled to 
biophysically detailed models of cell contraction. 


3.2.1 Fundamental Equations 


In this study we will primarily use the active stress approach, but for completeness 
we also present the equations arising from the active strain approach. This model 
can be derived as a direct extension of the single myocyte model in (14), using a 
similar approach as in to consider both the intra- and extracellular domains: 
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Fig. 3.2: Illustration of the intra- and extracellular domains for a single cell and its 
surroundings. 


a: V. Pj «0, b: Pi = Sg. c: F; = FPF?, in Qj, 

d: V-Pe=0, e: Pe= Se, in Qe, 

fi Ui = Ue, g: ni: Pi = ne: Pe, onT, (3.1) 
h: ne: Pe =0, on 0. T, 

i: u=0, on 0O,. p. 


Here, Q; and Qe are the intra- and extracellular domains, respectively, I is the inter- 
face between the domains, with the normal vector n; pointing out of the intracellular 
domain and n, out of the extracellular domain. Finally, OQ. r and 0€, p are the 
parts of the outer boundary Qe subject to traction- and displacement boundary con- 
ditions, respectively. See Figure 3.2 for a sketch of a typical computational domain, 
including a single cell and its immediate surroundings. Following (14), we here 
apply the active strain approach to incorporate active contraction of the myocyte, 
where the intracellular deformation gradient F; is decomposed as described above. 
The passive part is assumed to be hyper-elastic and derived from a strain energy 
function, see for instance (6) for details. A common choice for the active part is 
F? = diag((1 - y). (1 — y) |, (1 — y)*12), where y describes the fiber contraction 
and is a function of the cell activation state. For a more detailed introduction and 
discussion of active strain models, we refer to (1). 


The active stress model is the most widely used approach for modeling coupled 
active and passive mechanics on tissue level, and this is the approach we will employ 
in the subsequent numerical experiments. In the present context the active stress 
model involves a decomposition of the intracellular first Piola-Kirchhoff stress P; 
into a passive elastic part p. and an active part P7. The passive stress is derived 
from a strain energy function in the usual way, while the active stress is a function of 
the cell activation state. For the simplified model considered here we write the active 
stress as a function of time and the local fiber stretch 4, but the approach can easily 
be extended to include detailed biophysical models of the contractile mechanisms. 
The full active stress model may be written as 
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a: VP; 20, b: P; = Ft + P*(t,4), in Q;, 

c: V.P, 20, d: Po= $F, in Qe, 

e : Ui = ue, f: ni: Pi =ne: Pe, onT, (3.2) 
gine: P, =0 on 0Q¢.,7, 

h:u=0 on 0O, p. 


Both approaches treat the extracellular domain in the same way, as a passive hyper- 
elastic material governed by a strain energy function Ye. As given by (3.1) f-g and 
(3.2) g-h we assume continuity of stresses P;, Pe and displacements u;,ue across 
the cell membrane I’, implying that the membrane itself has no stiffness. The outer 
boundary ©, is assumed to be stress free, with Dirichlet conditions applied to parts 
of the boundary to avoid rigid body motion. Models for the active stress P? come 
in many forms, including simple phenomenological models as well as detailed bio- 
physical models of cell electro-mechanics [13). For the present study we apply 
a simple model where the active stress is derived from a (pseudo-) strain energy in 
the same way as the passive stress: 


a 9T 
GU 


(3.3) 


Here, V? is given by 
Tactive(t) 
a _ +active 
Y7 = 2 
where A = ||F'ei|| is the stretch in the so-called fiber direction (i.e. the main orienta- 
tion of the muscle cells), defined by the unit vector e1, and Tactive(t) is a prescribed 
function defining the active contractile force as a function of time. 


PL 


3.2.2 Specific Model Choices 


In this section we describe specific choices of the constitutive laws describing active 
and passive material properties in the models above, to arrive at a complete model 
that can be solved for the deformations and stresses. As noted above, we will in the 
following only consider the active stress model, given by (3.2). For the strain energy 
defining the passive stress-strain relationships we have applied a model from (19), 
which belongs to the family of models first presented by Guccione et al. (5). The 
same form of strain energy is used in the intra- and extracellular domains, but we 
allow the material parameters to be different. Both domains are modeled as nearly 
incompressible, with volume changes during deformations controlled by a penalty 
term. We have 


V; = Cj(e9i — 1) + «(Jn J — J +1) x EQ, (3.4) 
V, = Co(e2e — 1) + «(J1nJ —- J +1) xX E Qe, (3.5) 
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where Q;, Qe are functions depending on components of the Green-Lagrange strain 
tensor E — i(FTF - D 


2 2 2 2 2 
Q; -b; jE}, + bi (Ey, + E33 + Ex + Ex 
tbe, (Et, + E2, + E2, + E2)). (3.6) 


Furthermore C}, by. j, b; j, and bys j, for j = i,e are material parameters characteriz- 
ing the material's stiffness to the various strain modes, x is a penalty parameter that 
controls the volume changes, and J = det F. For a fully incompressible deformation 
we have J = 1, and in our nearly incompressible model we tune the parameter x to 
keep J = 1. 


In its most general form, the materials described by (3.4)-(3.6) are are transversely 
isotropic, which is a special case of orhtotropic materials. While an orthotropic 
material has different mechanical properties in three different directions, a trans- 
versely isotropic material is isotropic in planes normal to a characteristic direction. 
Passive cardiac tissue is known to behave as an orthortopic material (9), with the 
three directions dictated by the orientation and organization of the myocytes. How- 
ever, a transversely isotropic material is shown to be a good approximation, with 
material isotropy in planes normal to the fiber direction, the main orientation of 
the muscle cells. The details of the intra- and extracellular material behavior in our 
micro-structural model are less well-studied, and the degree of anisotropy has not 
been characterized. From the microstructure of the contractile apparatus occupying 
most of the intracellular space (see Figure 3.1) it is natural to assume anisotropic 
behavior, but the exact degree of aniostropy is not known. As a starting point, we set 
the intracellular material parameters to 


br; = 8, bte = 2, bee = 4. (3.7) 

For the extracellular space we assume isotropic material behaviour, setting 
bf e = bre = bfse = 1. (3.8) 
The bulk compressibility was set to x = 1000 kPa in both domains, while we explored 


different values of the scaling parameters C; and Ce, to be specified below. 


For the active stress model defined in (3.3) we have used a pre-computed transient 
tension Tactive(t) as shown in Figure 3.3. The curve was computed using the model 
of Rice et al. with default parameters, which outputs a normalized force. This 
value was then scaled such that the peak value reaches 2 kPa, giving a reasonable 
contractile stress for our application. 
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Fig. 3.3: Transient tension Tactive(t) over time (left), first computed in (13), then 
scaled to give values on a reasonable scale. In the intracellular domain the active 
tension is homogeneously set to this value; in the extracellular domain there is no 
such tension, implemented as being set to zero for all time steps. 


3.2.3 Numerical Methods 


The problem defined by (3.2) is solved with the displacement u as the primary 
unknown. To solve the system with the finite element method, it is convenient to 
formulate it as a single PDE defined over the entire domain Q = Q; U Qe. Such a 
formulation is not possible for the strong form of the PDEs, so we first need to derive 
the weak form of the equations. Starting with (3.2)a, we define a suitable vector 
function space V(O;) defined over the intracellular domain, multiply the equation 
with a test function v € V(Q;) and integrate by parts, to arrive at a weak formulation 


f P; - Vvdx — fo - Pi)v = 0. 
Qi r 


This equation is to be satisfied for all v € V(O;). Performing the same steps for the 
extracellular domain, and using the boundary condition (3.2)g on the outer boundary, 


we get 
; P.-Vvdx — fo. -Pew = 0. 
Qe T 


This equation should hold for all test functions v € V(Qe), where V(Q_) is a suitable 
space of functions defined over the domain Qe. Using similar arguments as in (15), 
we can define a function space V(Q) as the set of functions defined over Q that 
belong to both V(Q;) and V(Q_) and are continuous over I. With this definition, we 
may add the two weak forms above to obtain 


f p, Vids - f(n LX Pe-Wvdx~ | (ne - Pe = 0. 
Q; r Qe r 


Which is to be satisfied for all v € V(Q). Since n, = —n;, the surface integrals over 
I cancel because of (3.2)f. We can also use (3.2)e to define a single displacement 
field over Q, and we are left with the following weak form: Find u € V(Q) such that 


3 EMI Mechanics 35 


iy 
ely Z x 


Fig. 3.4: A: Volume element of one single cell; lines indicate cross section area. B: 
Cross section along longitudial direction of the cell. C: Volume element, 5 x 5 cells; 
lines indicate cross section area. D: Cross section along longitudial direction of the 
cells. 


J P. Vvdx = 0, (3.9) 
Q 


is satisfied for all v € V(Q). with P defined by (3.2)b and (3.2)d in the respective 
domains. 


3.3 Results 


We here present a number of numerical experiments to illustrate the general behav- 
ior of the models defined above. The code is implemented using FEniCS, and an 
archieved version of the code is available, see (16). 


For the simulations we used two different meshes; one representing a single cell 
and one representing a sheet of five by five cells, see Figure 3.4. Both meshes 
include subdomains defining the intra- and extracellular domains, separated by the 
cell membrane. To avoid rigid body motion, we keep a few points in the middle fixed. 
The rest of the boundary is kept unloaded to allow free contraction of the cells. 
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For each experiment we calculated the Green-Lagrange strain tensor E and the 
Cauchy stress tensor c', given by 


| JF’ PFT x €Q; 
— |\FI-'P.FT otherwise. 


On matrix form we can can write these out as 


Ey, En E43 7011 012 013 
E = | £E En Ez o = |021 022 023 
E31 E32 E33 031 032 033 


and for each of these we present plots for the first and the middle components, 
(E11, E2, 011, 022), which characterize strain and stress in the fiber and cross-fiber 
directions. 


Fig. 3.5: Tracking points, for which we evaluate functions of interest across various 
experiments. The points are uniformly distributed on a line from one corner to the 
middle, in the xy-direction, corresponding to the cross-section shown in Figure 3.4. 
Two of them are both located in the extracellular subdomain, and one should expect 
them to show different patterns than the three located in the intracellular subdomain. 


We first considered a single cell, and simulated contraction over a single cardiac 
cycle with homogeneous active force applied throughout the cell. For this simulation 
we chose parameter values Ce = C; = 0.5. The results are presented in Figure 3.6, 
where we observe that the deformation follows the expected pattern of a contraction 
in the longitudinal direction of the cell. Furthermore, in spite of the homogeneous 
applied active stress we see slight spatial variations in the deformation state, resulting 
from the discontinuity of active force across the cell membrane. 


Similar patterns are observed in the simulation of the sheet of 25 cells, shown in 
Figure 3.7. In this experiment the same active stress transient through the intracel- 
lular domain of all the cells, with the same material parameters. We still observe 
spatial variations in the deformation pattern — each cells is affected by mechanical 
deformation around them. 


We then considered two cases where we kept all parameters but one fixed, exploring 
the choices of material stiffness parameters C, and C;. The results are presnted in 
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Fig. 3.6: First and middle components of the Cauchy stress tensor c and Green- 
Lagrange strain E, for a single cell. The plots to the left shows values plotted over 
time, for the first 500 ms (out of 1000), following tracking points as shown in Figure 
3.5. The plots to the right shows values over a cross-section as shown in Figure 
3.4, as the active tension reaches it's peak value. The grey rectangle indicates initial 
configuration. 


Figures 3.8 and 3.9. These simulations were again performed on a mesh representing 
a single cell, with active force applied as described above. For the first experiment 
we kept C, fixed at 0.5, changing Cj; that is, we let the material stiffness in the extra- 
cellular domain remain the same while increasing the stress/strain scaling parameter 
in the intracellular domain. As C; increases the material becomes stiffer, and for the 
same active stress applied, one should expect less contraction. This can indeed be 
observed; both components of the Cauchy stress tensor (in magnitude) and the strain 
tensor decreases everywhere, and for the last three parameter choices there is almost 
no difference in deformation. On the other hand, we still apply an active stress in the 
intracellular domain, and we observe that the strain close to the membrane doesn't 
change much even if it changes everywhere else. 


For the next experiment we changed to keeping C; — 0.5 constant, while increasing 
Ce. We can observe higher Cauchy stress for the first component, and lower Cauchy 
stress for the second component, with increasing values of Ce. The strain decreases 
for both components. This is exactly as expected — in one end of the spectrum, having 
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Fig. 3.7: First and middle components of the Cauchy stress tensor c', and and Green- 
Lagrange strain E, for 5 x 5 cells. Values are plotted over the cross-section as shown 
in 3.4, as the active tension reaches it's peak value. The grey rectangle indicates 
initial configuration. 


Ce = 0.5, one would expect the extracellular subdomain to not affect the intracellular 
domain as it's rather flexible. For a given tension in the intracellular domain, it will 
just move along quite easily, while the overall behaviour in the whole domain is 
governed by the contraction inside the cell. As C; increases, the material is modeled 
as stiffer and hence constrain the movement more. For very high values the material 
is so stiff that it hardly moves, efficiently keeping the membrane close to fixed. 


3.4 Discussion 


We have presented a general framework for modeling cardiac mechanics on a sub- 
cellular scale, by extending a model of the type defined in to the extracellular 
domain. A series of preliminary numerical experiments demonstrate that the model 
behaves as expected, with the discontinuity across the cell membrane giving rise 
to spatially varying deformation fields even though both the active stress and other 
model parameters are spatially homogeneous over the intracellular domain. 


The main purpose of this work was to present the model framework and to illustrate 
the general behaviour of the model, while more detailed investigations and model 
extensions are left for future studies. A complete list of model limitations and poten- 
tial extensions would be too extensive to present here, but it is worth commenting on 
a few of the most obvious ones. First, the model derivation above included a number 
of simplifying assumptions on the mechanical properties of the cell membrane and 
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Fig. 3.8: First and middle components of the Cauchy stress tensor c and Green- 
Lagrange strain E, for a single cell, as we vary the parameter C;, which defines the 
stiffness of the material in the intracellular domain. Panel A shows spatial variation 
over a cross-section of the cell (see Figure 3.4), at peak. Panel B shows how the 
value, at peak, changes in given tracking points (see Figure 3.5). 
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Fig. 3.9: First and middle components of the Cauchy stress tensor c and Green- 
Lagrange strain E, for a single cell, as we vary the parameter Ce, which defines the 
stiffness of the material in the extracellular domain. Panel A shows spatial variation 
over a cross-section of the cell (see Figure 3.4), at peak. Panel B shows how the 
value, at peak, changes in given tracking points (see Figure 3.5). 
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the two domains. The continuity of stress across the cell membrane implies that 
the membrane itself has no stiffness, which is obviously incorrect, but it may be 
a reasonable assumption for many applications. The impact of different membrane 
mechanical properties should be explored further in a future study. Similarly, both 
the intra- and extracellular domains are assumed to be hyperelastic materials, which 
is probably a fairly crude approximation of the actual behaviour. In reality both of 
these domains are complex compositions of fluids and various embedded proteins 
structures, and the material behavior is most likely quite complex. Visco-elastic ma- 
terial models could potentially be a more accurate description than the hyper-elastic 
models applied here, but the required level of detail and material model complexity 
remains to be determined. Finally, we have here assumed that both domains are 
initially in a stress-free resting state, while experiments have shown that the extra- 
cellular matrix shrinks considerably when the myocytes are removed. Thus indicates 
that the resting state is actually an equilibrium state with non-zero stress in both 
domains, and accurately capturing the overall mechanical behaviour may require 
including this pre-stress in the model. 


In general, the level of detail and complexity of the model formulation will be dictated 
by the application. Some applications may require further development of the model 
along the lines suggested above, while for studies of a more qualitative nature the 
simplest version would be sufficient. One obvious application of the developed model 
framework, where a fairly simple model would probably give interesting results, is 
to study the impact of heterogeneities in calcium concentration and mechanical 
properties on the contractile properties of cells and tissue. 
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Chapter 4 


Operator Splitting and Finite Difference 
Schemes for Solving the EMI Model 


Karoline Horgmo Jeger!, Kristian Gregorius Hustad!:?, Xing Cait? and Aslak 
Tveito!-? 


Abstract We want to be able to perform accurate simulations of a large number of 
cardiac cells based on mathematical models where each individual cell is represented 
in the model. This implies that the computational mesh has to have a typical resolution 
of a few um leading to huge computational challenges. In this paper we use a certain 
operator splitting of the coupled equations and show that this leads to systems that can 
be solved in parallel. This opens up for the possibility of simulating large numbers 
of coupled cardiac cells. 


4.1 Introduction 


In recent publications we have shown that a cell-based model is useful 
for accurately representing the electrophysiology of excitable cells. Traditionally, 
excitable tissue is simulated based on homogenized models where the cells are not 
explicitly resolved, see e.g., (261[7). In the cell-based model, we explicitly represent 
both the extracellular space (E), the cell membrane (M) and the intracellular space 
(D, and it is therefore referred to as the EMI model. Similar approaches to modeling 
excitable tissue have been used by several authors; see e.g., 
BA. 


The EMI model is solved, numerically, using an operator splitting scheme which 
results in two steps; a non-linear system of ordinary differential equations (ODEs) 
to be solved in each computational node (i.e, degree of freedom) placed on the cell 
membrane, and a linear system of algebraic equations coupling the discrete Laplace 
equations of E and I with continuity requirements of the current over M. The spatial 


1 Simula Research Laboratory, Norway 
?Department of Informatics, University of Oslo, Norway 


The Author(s) 2021 44 
A. Tveito et al. (eds.), Modeling Excitable Tissue, Simula SpringerBriefs 
on Computing 7, https://doi.org/10.1007/978-3-030-61157-6 4 


4 Operator Splitting and FD Schemes 45 


resolution used in the discretization of the model is usually between 1 um and 4 um, 
thus only 1 mm? of tissue leads to more than 107 computational nodes. For an adult 
human cardiac cell, with a resolution of 2 um, the number of computational nodes 
per cell (including the associated extracellular space) is about 6000 (see (30), Table 
7). Thus, for a limited number of cells, the linear system coupling all the discrete 
Laplace equations is manageable. In fact, the system was solved using Matlab for up 
to 16,384 cells, with about 9.8x 107 computational nodes, see (30). 


However, not only the sheer size of the linear system is a challenge, also the properties 
of the linear system are unusual. In scientific computing, one of the most well-studied 
problems is solution of linear systems arising from the discretization of elliptic 
boundary value problems; see e.g., (5] 21} [8). Unfortunately, the EMI system does 
not naturally fall into the category of elliptic boundary value problems that can be 
solved using well-developed numerical machinery. It is therefore of importance to 
develop a splitting strategy for the EMI model that leads to sub-problems of the 
elliptic type. In (14), we showed that such a splitting can indeed be achieved. Here, 
we will review this convenient way of splitting the EMI model and show how to solve 
the system numerically using a finite difference method. Moreover, we will use the 
numerical scheme to assess the conduction properties in a small collection of cells 
where a sub-group of the cells are ischemic. Furthermore, we will present a parallel 
implementation of the splitting strategy, based on using open-source numerical 
libraries. This code is considerably faster than the existing Matlab code, and well 
suited for shared-memory parallel computers. 


4.2 The EMI Model 


We model the electrical properties of collections of cardiac cells using the EMI 
model introduced in [30). In Figure 4.1 we show the computational 
domains in the case of two coupled cells. Here, en and Q? denote the intracellular 
domains, and Qe denotes the extracellular space. The cell membranes are denoted 
by I and I5, respectively. The intercalated disc at the intersection between en and 
qs allowing for currents between the cells, is denoted by I'; ». With this notation at 
hand, the EMI model takes the following form: 


V. oiVuk -0 in ak, Ne ` Te Vue = -nk . a; Vut = p at Ik, 
V.o,Vu, 20 inQ,, vk = ec s - E) at I, 
ue —-0 at AQP, ib eu = wk at Lid 
Ne OeVue —0 at aan, nk . oiVuk = =n . ca; Vut mii at T. g 
uk — ue =v" atTy, wk = a Up Top) at T, ġo 


sk = F* at I. 
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Fig. 4.1: A: Two-dimensional version of the EMI model domain in the case of two 
connected cells. Here, the cells o and Q2, with cell membranes denoted by I'; and I5, 
respectively, are connected to each other by the intercalated disc, I’) 2, and surrounded 
by an extracellular space, denoted by Qe. B: Two-dimensional illustration of the 
geometry used for a single cell. The intracellular domain of each cell is composed of 
five subdomains Oo, Ow, Og, Qs, and On. The sizes of the subdomains are specified 
in Table 4.1. 


The model is stated for cell number k, and k denotes one of the six neighboring cells 
(in 3D: north, west, south, east, above, below). In the model, ue, uk ,and v^ = uk —Ue 
denote the extracellular, intracellular, and transmembrane potentials, respectively. 
Also, w* is the potential difference across the intercalated disc!, T, ;, and oj and Fe 
denote intracellular and extracellular conductivities, whereas C, and C, represent 
the specific capacitance of the membrane and the intercalated disc, respectively. 
Furthermore, ne, nk , and nk represent the outward pointing unit normal vectors of 


Qe, Q% and OF, respectively. A homogeneous Dirichlet boundary condition is applied 
at the outer extracellular boundary in the x-direction (0Q?), and a homogeneous 
Neumann boundary condition is applied at the boundary in the y- and z-directions 
(QN). The parameters used in the computations below are summarized in Table 
4.1. The properties of the cell membrane and the gap junctions are represented by F, 
lion and Igap. In the computations reported below, we use the Grandi et al. model(9), 
to model the dynamics of the membrane (F and Jion), and for the gap junctions we 
use the simple passive model is = wk /Rg. 


4.2.1 Operator Splitting Applied to the EMI Model 


As mentioned above, a key step in solving the EMI model is to split the equations into 
parts that can be solved using standard tools. In (14), we derived a splitting scheme 
that leads to two key numerical challenges: Non-linear systems of ODEs to be solved 


1 Note that w* is defined specifically for each cell. 


4 Operator Splitting and FD Schemes 


Parameter Value Parameter Value 

Size Qo 100umx18umxi8um Cm 1 uF/cm? 
Size Ow, Og 4 um x 10 um x 10 um C, 0.5 uF/cm? 
Size Ou, Os 10 um x 4 um x 10 um Ci 4 mS/cm 

Ax, Ay, Az 2 um oe 20 mS/cm 

At 0.02 ms Rg 0.0045 kOcm? 
AtopE 0.001 ms Mi, Ni 2 


47 


Table 4.1: Parameter values used in the simulations, based on (13). For parameters 


of the Grandi model, see (9). 


Algorithm 1: Summary of the splitting algorithm for the EMI model for 
connected cells. 


Initial conditions: v^, 5.9, w*-°, u? for all k. 
forn=1,...,N;: 


time step At from (skort! ynl) of 


k 1 k ok 
Vi = -T lion >S ) 


sk = F(v*, s^). 


n-l wk = wln-l 


e 3 
forj=1,..., Ni: 
Step 2: 

form= 1,..., Mg: 


Define üe = u 


For every k, find ak by solving 


V-o:Vak = in OF, 

ay + Bonk - oi Vat =y" +e at Ix, 
k -k | 1 >k, wk-wk-n-l E 
-ni -o; Vu; = RU + Cg A: at D, i. 


where denotes each of the neighboring cells of cell k. 


Update w^ = üF -= üt at T} g for all k and k. 
end 
Step 3: Find à, by solving 


V -Oe Vite = 0 in Qe, 
ü,-0 at oQD, 
Ne + Te Vile =0 at oa’, 


Ne ` Oe Vile = -nk . c; Vik at I. for all k. 


end 


Define u? = ite, un? = ak, w* = w for all k. 
Step 4: Define v^^ = ul" - uz at y for all k. 


i 


end 


Step 1: For all k, find s*-” and v* at the nodes of the membrane T% of cell k by solving a 
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at each computational node located at the cell membranes, and a series of elliptic 
equations; see Algorithm 1. All the differential equations involved in Algorithm 1 
are of classical type and can be solved using well-established numerical methods. In 
our present implementation, we apply a straightforward finite difference scheme (see 
e.g., for an elementary introduction to finite differences) for the elliptic equations 
and the Forward Euler method with a substepping time step Afopg for solving the 
non-linear ODEs modeling the membrane dynamics (see (30)). However, it is worth 
observing that elliptic equations can as well be solved using the finite element 
method, or a finite volume method, thus allowing for more flexible and adaptive 
meshes. 


4.3 Simulating the Effect of a Region of Ischemic Cells 


In order to demonstrate an application of Algorithm 1 above, we consider a collection 
of cells where a fraction of the cells are ischemic. This is known to perturb the electri- 
cal conduction and may lead to arrhythmias; see e.g., (33][281[19] 27). This problem 
has been carefully studied using homogenized models (mostly the monodomain 
model), but here we will show that the ischemic regions also have local effects when 
only very few cells are considered. In Figure 4.2, we consider a collection of cells 
organized in a two-dimensional mesh of 22x12 cells. The cells are modeled using 
the Grandi model with parameters as stated in (9). Within the domain, 8x6 of the 
center cells are ischemic in the sense that the extracellular potassium concentration 
surrounding these cells is increased from 5.4 mM to 10 mM. For the ischemic cells, 
we use the steady-state values of the state variables for the increased extracellular 
potassium concentration as initial conditions, and for the remaining cells, we use the 
steady-state values of the default Grandi model. In addition, we run the simulation 
for 5 ms before stimulation. 


In the simulation results we observe that the ischemic region slows down conductions 
and thus perturbs the wave in the intracellular potential moving from left to right. This 
is consistent with the result obtained in (6) (cf. Figure 5), where the monodomain 
model was used. Such perturbations are known to be arrhythmogenic and have been 
Observed several times in numerical experiments; see e.g., [15] [31 [6). Here, 
we observe that such perturbations can be initiated locally when only a few cells are 
subject to surroundings with elevated potassium concentration. 
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Fig. 4.2: Extracellular potential (left) and intracellular potential (right) at four differ- 
ent points in time in an EMI model simulation with an ischemic region in the center 
of the domain, marked by the purple rectangle. The parameter values used in the 
simulation are given in Table 4.1. 


4.4 A Scalable Implementation of the Splitting Scheme 


In expectation of future simulations of excitable issues that may involve a huge 
number of cells, we see the need of a scalable implementation of the new splitting 
scheme, so that it can run efficiently on parallel computers. One specific criterion is 
that the computation time should grow linearly with the number of cells involved. 
Additionally, the design goals of this new code should also include independence of 
proprietary software (such as Matlab) and plug-and-play of the different numerical 
components. This section will present a preliminary version of such a scalable 
implementation. 


50 Jeger et al. 


4.4.1 The Linear System for the Intracellular Potential 


The main benefit of the new the splitting scheme is that the intracellular Laplace 
equations (one per cell) are decoupled from the extracellular Laplace equation, as 
stated in Algorithm 1. If we assume a constant intracellular conductivity o; and 
that each cell is of the same shape and size, as shown in Figure 4.1, the matrices 
arising from a standard finite difference discretization of the intracellular Laplace 
equations for the individual cells will be mostly identical. There are only a small 
number of unique intracellular matrices, depending on whether there is a neighbor- 
ing cell connected to each of the intercalated discs. It is thus unnecessary to compute 
an intracellular matrix for each cell. Instead, the cells that have the same neighbor 
connectivity situation can share the same intracellular matrix. This not only reduces 
the memory usage of an implementation, but also improves data reuse in the caches 
of a computer. Moreover, since the number of computational nodes per intracellular 
domain is relatively small (each intracellular domain has about 5300 degrees of free- 
dom for the simulations used in this chapter), it is very efficient to use a direct solver 
each time an intracellular Laplace problem needs to be solved. Specifically, the LU 
factorization of each unique intracellular matrix A; can be pre-calculated, which ren- 
ders the solution of A; ak = bk per cell to be merely invoking the forward-backward 
substitution procedure. Parallelism of the computation mainly arises from the fact 
that the intracellular Laplace equations can be solved independently of each other, 
while limited parallelism also exists within each forward-backward substitution. 


4.4.2 The Linear System for the Extracellular Potential 


For the overall extracellular Laplace problem, which can be huge depending on the 
spatial resolution and the total number of cells, an iterative solver is more appropriate. 
Take for instance the case of 128 x 128 cells. The corresponding discrete extracellular 
Laplace equation has 107,202, 214 degrees of freedom. Independent of the spatial 
resolution and the number of cells involved, the extracellular matrix Ag arising from a 
standard finite difference discretization is symmetric and positive-definite (some care 
is needed to discretize the boundary conditions on the membranes). The resulting 
linear system Agüe = be is thus a perfect candidate for the conjugate gradients 
(CG) method with an algebraic multigrid (AMG) preconditioner. Under optimal 
conditions, an AMG preconditioner requires a constant number of iterations to reach 
convergence independent of the linear system size, although the number of grid 
levels inside the AMG preconditioner may increase with the system size. Parallelism 
readily exists in iterative solvers, with several software libraries providing parallel 
implementations of CG and AMG. 
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4.4.3 The Non-Linear ODE System for the Membrane Potential 


For solving the non-linear ODE system per computational node on the membranes, 
a straightforward and often very efficient numerical strategy is the Forward Euler 
method with a substepping time step Atopg. Since the non-linear ODE system on 
each membrane node is independent of the others, the ODE computation possesses 
the most ample parallelism. 


4.4.4 The Implementation 


The Python programming language has been chosen for the implementation, mostly 
because of its flexibility for interfacing with numerical software libraries written 
in performance-friendly languages such as C and C++. We have used the ctypes 
module from the standard Python library for this purpose. The choice of Python also 
simplified a partial translation from the existing MATLAB code developed in (14). 


We have chosen the SuperLU library for performing the LU factorization of 
the intracellular matrices and the subsequent forward-backward substitution, via the 
bindings that are provided by SciPy (32). For the extracellular Laplace equation, we 
have used the ViennaCL library for its implementation of CG and AMG. The CG 
iterations are by default configured to terminate when a tolerance of 10^? is reached. 
The AMG preconditioner has been configured to use the maximum independent 
set (MIS), see (4), as the coarsening algorithm and smoothed aggregation as the 
interpolation algorithm. For the ODE part, the Gotran automated code generator 
has been used to translate the Grandi cell model into C code, callable from the Python 
side. 


4.4.5 Parallelization 


The Python implementation currently relies on the adopted numerical libraries (Su- 
perLU and ViennaCL) for an implicit parallelization of the PDE computation through 
multi-threading. This form of parallelization suits for shared-memory parallel com- 
puters, such as laptops or servers that use multicore CPUs. Multi-threading of the 
ODE computation is also enabled by inserting OpenMP compiler directives into the 
C code that is generated automatically by Gotran. The advantage of this implicit par- 
allelization is that the user does not have to care about parallelization-specific coding. 
The downside is that all the computations have to run on a shared-memory system. It 
is possible to achieve the more general parallelization that targets distributed-memory 
parallel computers, which will be a task for future work. 
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4.4.6 Performance Results 


The simulations in this section were run on a dual-socket server with two 32-core 
AMD EPYC 7601 CPUs, each with 8-channel memory operating at 2666 MT/s. 
The number of OpenMP threads was set by default to the number of logical cores, 
equaling 128. Moreover, the environment variable OMP_PROC_BIND=TRUE was set 
to prevent the threads from migrating between the cores (which typically leads to 
unnecessary performance loss). 


Table 4.2 shows the average solution time per time step for the 10 first time steps, 
where all the parameters are as prescribed in Table 4.1. The number of cardiac cells 
is doubled in the x and y directions for each row, and we observe that the time per 
cell remains fairly constant, indicating that the time to solution is a linear function 
of the number of cardiac cells simulated. 


Cells time usage for all cells (s) time per cell (ms) 
E M I total E M I total 
4x4 0.38 0.03 0.09 0.50 24.0 22 53 31.5 
8x8 1.54 0.11 0.34 1.99 241 18 52 312 
16x 16 2,27 0.45 1.20 3.92 8&8 17 47 153 
32 x 32 8.91 1.72 4.98 15.61 87 17 49 15.2 
64 x 64 30.46 6.73 19.15 56.33 74 16 47 13.8 


128x128 123.73 30.60 72.83 227.16 76 19 44 139 


Table 4.2: Average solution time per time step for the E, M and I domains. 


4.5 Software 


The Matlab code used to compute the solutions shown in Figure 4.2 and the Python 
code discussed in Section 4.4 can be found at https: //github.com/KGHustad/ 
emi-book-2020-splitting-code| An archived version is also available. 


4.6 Conclusion 


In this chapter we have presented a numerical scheme for solving the EMI equations 
using operator splitting. The scheme allows for parallel solution of individual cells 
combined with a global solution of the equation modeling the extracellular potential. 
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The latter is well suited for using optimal linear solvers such as AMG. The overall 
code scales linearly with the number cells and thus allows for simulation of a large 
number of cells. It remains to be seen how well this will work for very large numbers 
of cells; this is subject for further work. 
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Chapter 5 


Solving the EMI Equations using Finite Element 
Methods 


Miroslav Kuchta!, Kent-André Mardal'? and Marie E. Rognes! 


Abstract This chapter discusses 2 x 2 symmetric variational formulations and as- 
sociated finite element methods for the EMI equations. We demonstrate that the 
presented methods converge at expected rates, and compare the approaches in terms 
of approximation of the transmembrane potential. Overall, the choice of which for- 
mulation to employ for solving EMI models becomes largely a matter of desired 
accuracy and available computational resources. 


5.1 Introduction 


In this chapter, we present different weak formulations and corresponding finite 
element methods for solving the EMI equations as presented in (7, Chapter 1) over 
a physiological cell Q; and its membrane I surrounded by an extracellular space Qe 
and a time interval (0, T] for some time T > 0. This coupled, time-dependent, and 
typically nonlinear system of equations can be targeted numerically by an operator 
splitting scheme, see e.g (8| Chapter 4). Such an approach, combined with for instance 
an implicit Euler discretization in time, gives the following stationary and linear, but 
still coupled system of equations to be solved at each time-step: find the potentials 
Ue = ue(x), uj = uj(x) (and current Im = Im(x)) such that 
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-V-. (ae, Vue) 20 in O,, (5.1a) 
-V -(a;Vu;) 20 in Q; (5.1b) 
OeVue - Ne = —0; VU; - Ni = Im onT, (5.1c) 
Uj —Ue = V onT, (5.1d) 

y - CAL, = f on I, (5.1e) 


where At > O denotes a time step size, and ne (resp. n;) denotes the outward 
pointing normal on I when viewed from Qe (resp. Q;). In our (implicit Euler) time 
discretization context, the known right-hand side f of (5.1e) combines the previous 
transmembrane potential solution, vo, and the evaluation of the ionic current, Jion, 
into f = vo — Cr At hone 


We assume that the potential is grounded on part of the external boundary T? and 
that the remaining external boundary I is insulated. These assumptions give the 
boundary conditions: 


Ue =0 on IP, (5.2a) 
D, Vu, ne =0 on TY. (5.2b) 


This geometrical setting is illustrated in Figure 5.1. 


Fig. 5.1: (Left) Illustration of the geometric setting for the single cell EMI problem. 
(Right) Sample meshes for our benchmark problem (5.17). The boundary facets of 
the intracellular mesh Q;,„ form the membrane mesh T}. 


Remark 5.1 We remark that a single cell model is here considered for simplicity. In- 
deed the formulations to be studied below can be similarly derived for the intercalated 
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model (collections of connected cells). Formulations for a number of disconnected 
cells are then practically identical to the case considered here. 


Remark 5.2 If (5.2) is considered without any Dirichlet boundary data, i.e. |T? | = 0, 
then only the transmembrane potential is fixed and the intracellular and extracellular 
potentials are determined only up to a single, common constant. 


The EMI equations (5.1) set a rich scene for numerical exploration and can be solved 
in a multitude of ways. In this chapter, we will derive 2x2 different weak formulations 
(each defining a finite element method) of this system. The two first formulations 
(in Section 5.2) compute the intracellular and extracellular potentials as the main 
unknowns. These are referred to as primal formulations. The latter two formulations 
(in Section 5.3) additionally introduce the current densities J; = —c; Vu; and Je = 
—Ge Vue as independent unknowns. These are referred to as mixed formulations. We 
compare finite element discretizations of the primal and mixed formulations with 
respect to the approximation of the transmembrane potential v in Section 5.4. This 
choice is motivated by the observation that v is closely coupled to the membrane 
dynamics as discussed in Chapter 1. 


5.1.1 Preliminaries: Function Spaces and Norms 


The EMI equations (5.1) define a multi-dimensional! PDE system coupling unknown 
fields defined over cellular domains and fields defined over the cell membrane, which 
can be viewed as a lower-dimensional manifold. Identifying the right function spaces 
for the different unknown fields is key to defining well-posed weak formulations of 
these equations. We here present suitable Sobolev spaces for this setting; the reader 
is referred to e.g. for more material and careful formalizations. 


Let Q be a bounded, polygonal domain in R for d = 2,3. We denote the space of 
square-integrable functions over Q by L?(Q), and let H!(Q) be the Sobolev space 
of functions in L^(Q) with weak derivatives in L*(Q). The space H (div, Q) contains 
vector-valued functions v : Q — R? such that v € L?(Q) and V - v € L?(Q). In 
general, when clear from the context, the domain will be omitted from the notation. 


The L?-inner product and norm for u,v € L?(Q) is written as 


2 2 
(na = f wae, Iva = f» dx. 
Q Q 


Similarly, we define the H!-norm as Ivl? o a llvil ot IVI o for v € H!(Q), and 
the H (div)-norm as lees - llvllo o + [V vlo 


1 PDEs coupling fields over domains of different topological dimensions are often referred to as 
mixed-dimensional PDEs. To avoid the confusion-inducing term mixed-dimensional mixed in the 
subsequent sections, we instead use the term multi-dimensional in this chapter. 
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For IT € ðQ, we define the constrained spaces Hy(Q) E {v e H'(Q)|v =Oon r}, 
and Hr(div, Q) = (v € H(div, Q) | v. n = 0 on T} where n is the (outward pointing) 
normal vector of T. Finally the spaces H'/?(T) and H~!/?(L) are the trace spaces of 
H! and H(div) respectively (6] Ch. 1., 2.). Here, the spaces will be considered with 
the norm defined in terms of fractional powers of the Helmholtz operator, see e.g. 
(4), i.e. 

lll? = (u, (—A  D'ior. u € C" D). 


We remark that in the following experiments the fractional norm is evaluated using 
the eigenvalue decomposition of —A + J as detailed in (11). 


5.2 Primal Formulations 


We present two primal formulations of the stationary EMI system (5.1) with the 
boundary conditions given by (5.2): one single-dimensional formulation and one 
multi-dimensional formulation. The difference in the intra- and extracellular potential 
across the cell membrane [I sets up a potential jump, the transmembrane potential 
v, c.f. (5.1d). Due to this jump, one cannot define a global, differentiable potential 
u € H'(Q; U Qe) such that u|o, = ui and u|o, = ue. Instead, we seek u; € H'(Q;) 
and ue € H!(Q,) separately. In the single-dimensional formulation, these are the 
only unknown fields, while in the multi-dimensional formulation, we keep Im as an 
additional unknown. 


5.2.1 Single-Dimensional Primal Formulation 


Define the function spaces 
V 2H(Q), Ve = Hep (Qe). (5.3) 


To derive a weak formulation of (5.1), multiply (5.1a) by a test function ve € Vo, 
(5.1b) by another test function v; € V;, and integrate the divergence by parts. This 
yields the variational formulation: find ue € Ve, u; € V; satisfying 


f TeVe : Vve dx — ] 7» : NeVe ds = 0, (5.4a) 
[o Tr 


f oa, Vu; : Vv; dx + ] coo : nj) vids = 0. (5.4b) 
Q; r 


i 


for all ve € Ve, vj € V;. In the bracketed term of (5.4b), we recognize the membrane 
current density 7,, as defined by (5.1c), and similarly, the interface contribution in 
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the corresponding extracellular equation (5.4a) hides —J,,. Combining (5.1e) and 
(5.1d), we find that 

Im = C(t)! (ui — ue) - f). (5.5) 
After substituting (5.5) into (5.4), the single-dimensional primal weak form of (5.1) 
reads: find u; € V; and ue € V, such that 


| CeVue ' Vve dx + S 6e. ds— i Cin(At) uive ds = 
Qe T T 


7 (Any! f fveds, 

i (5.6) 
1 ciVui : Vvi dx f Chry us ds- f Cn)" iv ds = 
Qi T T 


i 


Cm (Ar) | i f vids, 
r 


for all ve € Ve and v; € V;. 


We remark that (5.6) can be viewed as a coupling of two Poisson problems with a 
Robin boundary condition on T. The well-posedness of the problem is then discussed 
in Chapter 6). Finally, the transmembrane potential can here be computed from 
its definition (5.1d) as a difference of the computed potentials. 


5.2.2 Multi-Dimensional Primal Formulation 


An alternative formulation can be derived by keeping /,, as a separate unknown 
field. Since T is of a different (lower) dimension than Q;, Qe; and as Im : I —^ R 
while u; : Q; — R, ue : Qe  R, we will refer to this as a multi-dimensional 
primal formulation. Observe that (5.4) now yields two equations for three unknowns 
uj € Vi, Ue € Ve, and Im € Q: 


1 OeVue : Vve dx — finve ds =0, Vve E Ve, 
Qe r 
f o;Vu; : Vv; dx + f» ds=0, Vv € Vj. 
Qi r 
However, the missing equation can be obtained from (5.5). Let 
Q-H'(T, Qt =H P0). (5.7) 
We remind the reader that if I'is a co-dimensional 1 subset of © then trace operations 


from Q to T, Tu = ulp, u € C(Q) and T,T = T|r : n, tT € C(Q), formally have 
the following mapping properties T : H!(Q) — HT) and T, : H(div,Q) > 
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H-'/?(L). Hence, let jm be aa test function from Q*. We shall then enforce that (5.5) 
holds in the weak sense: 


fo = Ue)im ds zi [aca tn ds = n ds, Vim € Q'. 
Tr T T 


In turn, the multi-dimensional primal formulation of (5.1) reads: find u; € Vj, 
Ue € Ve, Im € Q* such that 


f Oe Vue : Vve dx — f ds = 0, 
QA Tr 


| o; Vu; : Vv; dx + f vilg ds = 0, (5.8) 
Q; T 


[nts f winds- f eis = f fin ds, 
T T T T 


for all v; € Vi, ve € Ve and jm € Q*. We remark that (5.8) is closely related to 
the Babuska problem for enforcing boundary conditions by Lagrange multipliers 
and the mortar finite element method, see e.g. (13). With regards to evaluation 
of the transmembrane potential, we note that v can be post-computed in several 
ways: from (5.1d) (as for the single-dimensional primal formulation (5.6)) or from 
Im and (5.1e). 


5.3 Mixed Formulations 


We now turn to consider mixed formulations of the EMI system (5.1). Let us 
(re)introduce the current densities 


Ji = —0; Vti, Je = —Oe Vue (5.9) 


and the global field J on Q = Q; UQ, such that J|o, = J; and J|o, = Je. In general, 
we use the convention that for a scalar or vector field u defined on Q, the restriction 
on Q; and Qe is denoted by u; and ue, respectively. 


With these definitions, (5.1a)-(5.1c) become: find the current densities J;, Je (or J) 
and the potentials u;, ue (or u) satisfying 


oz! Je + Vu, 20 on Qe, (5.102) 
o Jit Vu; 2 0 on Q;, (5.10b) 
-V-J=0 in Q, (5.10c) 


Je: ne + Ji: nj =0 on I. (5.10d) 
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We refer to (5.10) together with (5.1d)-(5.1e) as the mixed EMI system with bound- 
ary conditions given by (5.2). Weak formulations of the mixed form can enjoy 
improved conservation properties and stability properties (3). In particular, approxi- 
mations of J may be computed such that they are exactly divergence free, cf. (5.10c). 


The continuity condition (5.10d) ensures that the normal component of J is contin- 
uous on I’. We remark that v € H (div, Q) implies continuity of v - n on I. Moreover, 
we observe that (5.10c) involves only divergence of the field J. It is therefore suffi- 
cient to seek J in S = Hy (div, Q). Note that in contrast to the primal formulation, 
here the Neumann boundary condition (5.2b) is enforced as an essential condition; 
that is, it is included in the construction of the function space S. 


5.3.1 Single-Dimensional Mixed Formulation 


Let 
S= \J € Hyy (div,Q); J -n € Lo}, v = LX). (5.11) 


To derive a weak form of the mixed EMI system, consider a test function T € S. 
Taking the dot product of (5.102), (5.10b) with T;, Te, integrating and applying 
integration by parts then yields 


[i oe tedx- f iS rude furo neds =~ f UeTe ` Ne ds, 
Qe Qe r r? 
| Jonas- f nV iridye [uri mids =0. 
Q; Q; T 


Observe that by continuity of the normal component of the test function (T; = Te:n 
on I), and the identity n, = —n;, the integrals on I can be added, resulting in 
no — ue)t - nj. Moreover, using (5.5), the membrane term can be rewritten as 
qe (Ch! AtJ : ni + f) t: ni. In turn, we arrive at the variational problem: find J € S, 
u € V such that 


fI rass [etaim nds f nae - f fenis 
Q r Q T 


- f a7- Jax =0, 
" (5.12) 


for all r € S and q € V, with c defined naturally as o|o, = c; and likewise for Qe. 
Note that due to the extra trace regularity of the trial/test space S all the terms in 
(5.12), and in particular the interface term Jr C AtJ - njt - nj ds, are well defined. 
Without the extra regularity, i.e. if $ = Arn (div, Q), this would not be the case. 


We remark that (5.12) is a I-perturbed mixed formulation of the Poisson problem 
(see e.g. [9) for more details on mixed formulations of the Poisson problem). 
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Considering the task of approximating the transmembrane potential, we observe that 
v can be computed in two ways, as for the multi-dimensional primal formulation. 
Indeed, in addition to the identity v = u; — ue, cf. (5.1d), equation (5.1e) can be used 
since Zm = J : nj is readily available. 


5.3.2 Multi-Dimensional Mixed Formulation 


As for the primal formulations, the multi-dimensional mixed formulation is obtained 
by keeping the interface term as an explicit unknown field. Let 


S = Hyn (div, Q), V =Q), W =H’). (5.13) 


To complete the formulation, the equation to be enforced weakly by test functions 
w € W is the membrane dynamics condition (5.1e) written in the form 


J -ni — Cm At) ly 2 -C,(Aty ! f onr. 


The final multi-dimensional mixed weak formulation then reads: Find the current 
densities J € S, potentials u € V, and transmembrane potential v € W such that 


fracas [nan fn mas eo, 
Q Q r 


- [av sax =0, (5.14) 
Q 


[rimas - f esty as - (Any! f fas. 
r r r 


for all r € S, q € V and w € W. Note that (5.14) is defined on the standard H (div) 
space, cf. (5.12), as the formulation no longer contains the troublesome interface term 
n C, At J - nit nj ds. With regards to the approximation of v in formulation (5.14), 
observe that no post-processing is required to obtain this quantity. This is contrast 
to the previous three formulations. We remark that (5.14) is closely connected to the 
Babuška problem for the mixed Poisson equation (2). 


5.4 Finite Element Spaces and Methods 


To solve the primal and mixed, single- and multi-dimensional weak formulations nu- 
merically, we approximate the continuous function spaces by discrete finite element 
spaces. Each choice of formulation and finite element space defines a finite element 
method for solving the EMI system. 
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Now, let €, be a mesh of the domain €) = Q; U Qe with characteristic mesh size A, 
which conforms to T in the sense that no element of Q, has its interior intersected 
by I. The meshes Qen and Q; n of the extracellular and intracellular domains are 
formed as non-overlapping subsets of the cells of QO;. As a consequence, the mesh 
Ij, of I is formed by the facets of elements of Oy, cf. Figure 5.1. We remark that the 
single-dimensional primal formulation allows for independent discretizations of Qj, 
Qe as well as T. 


The choice of the finite element spaces plays a crucial role for the stability of the 
different discrete formulations. In particular, for the saddle-point systems, the spaces 
must be compatible in the sense of Babuska-Brezzi and satisfy discrete inf-sup 
conditions, see e.g. (3). For the primal formulations (5.4) and (5.8), we seek discrete 
unknowns and test functions in 


Vin = P\(Qi,n) C Vi, Van = Pi(Qe,n) C Ve, Qn = Pih) C Q, (5.15) 


where P, denotes the space of continuous piecewise linears (defined relative to the 
relevant mesh). With these spaces, we expect linear convergence with the mesh size 
h for all variables in H!-norms and quadratic convergence in the L?-norm. 


For the mixed formulations (5.10) and (5.12), we seek discrete unknowns and test 
functions in 


Sn = RT)(Qn) C S, Va = Po(Qn) CV, Wy = Polha). (5.16) 


Here RT, denotes the lowest order Raviart-Thomas finite element spaces and Po 
denotes the space of piecewise constants defined relative to the relevant mesh. These 
spaces satisfy the relevant stability conditions, and we expect linear convergence of 
all unknown fields in their respective natural norms. 


5.5 Numerical Results 
5.5.1 Comparison of Convergence between Formulations 


In order to compare the properties of the different formulations, and in particular 
their numerical stability and accuracy, we consider a manufactured solution test case 
with a smooth analytical solution. We define Q = [0, 1]? and Q; = [0.25,0.75]? with 
(| = 0. For simplicity, let c; 21,0, = 2, Cm = 1, At € {1, 1077) and consider 
the exact solution 


Ue = sin(z(x + y)), 
(5.17) 
uj = TE + cos (x(x -Da- p) cos (xo -io-i I 
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which corresponds to (At dependent) right hand sides f = uj — ue — Atl. Note, 
that with (5.17) both v + O and J, + 0. We discretize the domain by a uniform 
mesh by dividing the unit square into n x n squares and dividing each square by 
the (left) diagonal into isosceles triangles of size h, cf. Figure 5.1. To compare the 
dimensionality of the different formulations, Table 5.1 lists the dimensions of the 
four different finite element pairings over these meshes. 


h |Ve.n] IVi al lQn] ISa] IVa | IW | 


4.42E-02 864 289 64 3136 2048 64 
2.21E-02 3264 1089 128 12416 8192 128 
1.10E-02| 12672 4225 256 49408 32768 256 
5.52E-03| 49920 16641 512 197120 131072 512 
2.76E-03| 198144 66049 1024 787456 524288 1024 
1.38E-03| 789504 263169 2048 3147776 2097152 2048 
6.93E-04|3151872 1050625 4096 12587008 8388608 4096 


Table 5.1: Dimensions of the different finite element spaces for uniform refinements 
of the unit square. The first row corresponds to a mesh of Q with n = 16, i.e. having 
2 x 16 x 16 cells. 
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Aa 319 Folly = velki 3 
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[ = urfi: 7 102 L see ee J 
P " | 102 F E ou | 
10 E 2A. Au - unlo: 2 E E P Wt d | 
^ [u — ui]: 1 punc 
10900 7 A |lu — unllo: 2 E 5 eee | 
E i o Zn ae Ina i: 1 j 103 E pi i l E 
1025 10? 10-55 1025 10? 10715 
h h 


Fig. 5.2: Convergence properties of (left) primal formulations (5.6)—(5.8) and (right) 
mixed formulations (5.12)-(5.14). The EMI system (5.1) is solved with the exact 
solution given by (5.17) and At = 1. Filled symbols correspond to single-dimensional 
formulations. The number associated with each line indicates the convergence rate 
obtained from a least squares fit of the corresponding data. 
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On a series of meshes and for a set of different timesteps, we compute ? the different 
approximations of all four finite element methods. We then evaluate the approxi- 
mation error by evaluating the difference between (higher order interpolants of) the 
exact solution and the approximations in appropriate norms. Figure 5.2 shows the 
errors of the different formulations (with At = 1). In the primal formulations, u 
refers to the global potential, i.e. u; = ulo;, ue = u|o,. As the formulations seek 
the approximations u; € H'(Q,), ue € H!(Q,) the error is considered in the natural 
(broken) norm ||u||j = (||u; In +||Ue Ies. We observe that all the quantities converge 
linearly in their respective natural norms, as expected. In particular, the errors in the 
current density /,, in (5.8) and the transmembrane potential v in (5.14) are reported 
in the fractional norms H~!/? and H!/?, respectively. The former is computed by 
first interpolating the error into the space of continuous piecewise cubic polynomials 
on Ip while for v — v; the P; element is used for error interpolation. Without this 
higher-order approximation of the error, i.e. if the error is computed in the same 
space as the discrete solution, we observe quadratic convergence. 


Finally, we note that the primal formulations yield identical approximations of u 
cf. Figure 5.2 (left). Similarly, the mixed formulations give identical approximations 
of (u, J) cf. Figure 5.2 (right). Considering for comparison the error in the potential 
in the L?-norm, it can be seen that the primal formulations are more accurate 
than the mixed formulations. The same experiments for At = 107^ give similar 
approximation results. However, it is not true that these conclusions hold in the limit 
of At approaching 0, see e.g. Chapter 6. 


5.5.2 Post-Processing the Transmembrane Potential 


With the exception of the multi-dimensional mixed formulation (5.14), the trans- 
membrane potential v in the remaining EMI formulations is computed by post- 
processing. In (5.6) the approximation v;, can be obtained by interpolating the dif- 
ference ui, n — Ue,h onto e.g. the space of continuous piecewise linear functions over 
Ij. This procedure can, of course, be used in the other formulations as well. How- 
ever (5.8) and (5.12) also offer an alternative approach. In the multi-dimensional 
primal formulation, the discrete membrane current density, Zm,» is computed in the 
space Pı(T}) of continuous piecewise linear functions on Tp. In turn, vp can be 
computed (in the same space) as a projection of C, Ar I, ;, + f. The same formula 
can be applied in the single-dimensional mixed formulation since the current density 
can be evaluated as J, - n;. We recall that in (5.12) the natural space for vp is the 
space of (discontinuous) piecewise constant functions on lp however. 


Convergence of the transmembrane potential obtained by the different formulations 
and post-processing strategies is shown in Figure 5.3 for the same test case as previ- 


2 The code used to produce results in this chapter is available athttps://github.com/MiroK/ 
and archived at (12). 
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ously. We observe that the primal formulations yield quadratic convergence and that 
the single-dimensional primal (5.6), and multi-dimensional primal formulation (5.8) 
are practically identical. The discrete potentials obtained by solving the mixed for- 
mulations converge linearly with the projection method in the single-dimensional 
mixed formulation yielding the most accurate vy. In particular, the approximation 
is better than that of the multi-dimensional mixed formulation for this test case. 
Computing the potential in the single-dimensional mixed formulation (5.12) by in- 
terpolating uj; — ue, leads to poorer approximation than for the multi-dimensional 
mixed formulation. By comparing the results for two different time steps, we observe 
that the rates do not change considerably if At is modified. 
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Fig. 5.3: Approximation of the transmembrane potential by the different EMI for- 
mulations. (Left) At = 1, (right) At = 107^. Postprocessing by projection (using the 
current density) is indicated by Zm. In multi-dimensional mixed formulation vp is 
obtained by solving (5.14). Interpolation of ui h — ue, is used in other formulations. 
The final number indicates the convergence rate. 


5.6 Conclusions and Outlook 


All four finite element formulations provide a converging approximation to the sta- 
tionary problem (5.1). This system is a key building block in any operator splitting 
algorithm for the time-dependent EMI equations (1.30). The formulations provide 
solutions which differ by accuracy as well as computational cost, cf. Figures 5.2-5.3 
and Table 5.1. The formulations also differ in robustness of their approximation 
properties with respect to At. This issue is however beyond the scope of this chap- 
ter, and the interested reader is referred to the discussion in Chapter 6. In terms 
of coupling with the membrane dynamics the single/multi-dimensional primal and 
single-dimensional mixed formulations require post-processing. However, all ap- 
proaches discussed in Section 5.5.2 are easy to implement. Therefore, the choice of 
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which formulation to employ in solving the EMI model is largely a matter of desired 
accuracy and available computational resources. 


Open Access This chapter is licensed under the terms of the Creative Com- 
mons Attribution 4.0 International License (http://creativecommons.org/licenses/ 
by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in 
any medium or format, as long as you give appropriate credit to the original aut- 
hor(s) and the source, provide a link to the Creative Commons license and indicate 
if changes were made. 

The images or other third party material in this chapter are included in the chap- 
ter’s Creative Commons license, unless indicated otherwise in a credit line to the 
material. If material is not included in the chapter’s Creative Commons license and 
your intended use is not permitted by statutory regulation or exceeds the permitted 
use, you will need to obtain permission directly from the copyright holder. 
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Chapter 6 
Iterative Solvers for EMI Models 


Miroslav Kuchta! and Kent-André Mardal!-? 


Abstract This chapter deals with iterative solution algorithms for the four EMI 
formulations derived in Chapter 5). Order optimal monolithic solvers robust with 
respect to material parameters, the number of degrees of freedom of discretization 
as well as the time-stepping parameter are presented and compared in terms of 
computational cost. Domain decomposition solver for the single-dimensional primal 
formulation is discussed. 


6.1 Introduction 


Spatial discretization of EMI models describing a few cells with a complex/realistic 
geometry or a large collection of cells leads to linear systems with considerable num- 
ber of unknowns. In our largest simulations, we may have linear systems involving 
billions of unknowns at millions of time steps. It is the purpose of this chapter to 
discuss how such systems can be solved efficiently with available algorithms. 


Let us denote the system size as N and let h be a typical grid/mesh size. The complex- 
ity of a solution algorithm can then be analyzed in terms of how the computational 
time grows with N. For instance, a naive Gaussian elimination would perhaps scale 
as O(N?) and for linear system involving billions of unknowns such an approach 
is therefore not feasible during a life-time. Here, we shall aim for algorithms that 
are order optimal, i.e. their solution time scales linearly, O(N), with respect to the 
number of unknowns. Clearly, O(N) is order optimal in the sense that this is also the 
complexity of writing the results to file. In general, direct solvers, like for instance 
LU, do not scale linearly in N. Here, we will use Krylov solvers like the Conjugate 
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Gradient (CG) or the Minimal Residual (MinRes) methods. Both methods are known 
to provide efficient computations when combined with proper preconditioning tech- 
niques. Preconditioners based on multigrid and/or domain decomposition methods 
have been shown to be order-optimal for a variety of problems, but the theory for 
the EMI problem is currently limited. Efficient solvers for the EMI-problem are 
discussed in the following. 


We say that a solver for the transient EMI model (1.30) is order optimal if the solution 
time grows linearly in At~}, where At is the time step. Considering the stationary 
problem (5.1), which is to be solved at every step of the temporal loop, it is clear that 
order optimality of the transient solver requires that the solver of the linear system 
does not degenerate for small (or large) At, i.e. that it is robust with respect to the 
time stepping. 


In this chapter we discuss solution algorithms for the linear systems due to the finite 
element discretization of (5.1) which are robust in h and well as At. Two types 
of approaches will be considered. Monolithic approaches, where all the unknowns 
are solved for at once, are the subject of Section 6.2. Section 6.3 then concerns a 
domain decomposition approach where one iterates between the intra/extra-cellular 
subproblems. The solvers will be compared in terms of robustness and cost, however, 
only serial performance will be addressed. We remark that parallel scalable solvers 
suitable for the mixed formulations of the EMI models (5.12) (5.14) are the balancing 
domain decomposition methods, see e.g. (24][25). For the elliptic single-dimensional 
primal formulation (5.6), the FETI domain decomposition methods, e.g. 
could be used. Moreover, to simplify the exposition the focus shall be on a single cell 
model. We remark that all the algorithms presented further can be generalized in a 
rather straightforward manner to models with multiple disconnected cells. However, 
construction of robust monolithic solvers for multi-dimensional formulations for 
collections of connected cells is out of the scope of this manuscript. 


6.2 Monolithic Solvers 


Let 71,x, = Ln be a linear system due to discretization of the continuous problem 
Ax = Lin W’, where W’ is a dual space to some Hilbert space W and A: W > W”. 
Note that in case of (5.1) the continuous operator A depends on material parameters 
Ci, Te, Cm as well as the time step size At. Here the focus is placed on the latter 
dependence. 


The monolithic solvers considered here are preconditioned Krylov methods where 
the preconditioner 8 : W’ — W shall be constructed such that the number of 
iterations required for solving problems 8, Anxn = Bn Ln is bounded in h and At. A 
constructive framework for establishing 8 is operator preconditioning (19), where the 
structure of the preconditioner reflects mapping properties of A as an isomorphism 
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between suitably chosen Hilbert space and its dual space. In particular, if A : W — 
W' is an isomorphism then a Riesz map with respect to the inner product of W is a 
suitable preconditioner. As explained in (19), multilevel or domain decomposition 
algorithms are equivalent to Riesz maps for standard elliptic problems. Furthermore, 
these schemes also work in a fractional setting (1); a property which will be exploited 
in this chapter. We shall illustrate the framework briefly below using the single- 
dimensional formulations (5.6), (5.12) as an example. 


In the following we follow the notation of Chapter 5. In particular, ||-||o,o denotes the 
L? norm of a scalar or vector field over Q, while ||-||1,0, ll-llaiv,o are the H!(Q) and 
H (div, Q) norms respectively. For a Hilbert space W other than L?, H!, H (div) the 
norm is denoted as ||-||w while (-,-)w denotes the inner product on W. Moreover, we 
let (,-)o be the L? inner product. If the domain is clear the subscript will be omitted. 
By (.,-) we shall also denote a duality pairing between a Hilbert space and its dual. 
The dual of an operator A with respect to the L? inner product is then denoted as A'. 
Finally, let us introduce scaled, sum and intersection spaces, which will be required 
for well-posedness of some of the formulations. If W, Q are Hilbert spaces and a an 
arbitrary positive real number, the scaled space aW, the sum space W + Q and the 
intersection space W n Q are Hilbert spaces with norms ||x|law = all - Ilw, 


= 1 2 2 = 2 2 
lxlhwo = „int, Viele, Melo and [hxthwoo = Vli, + Ixl. 


weW.qcQ 


6.2.1 Single-Dimensional Primal Solvers 


In order to simplify the analysis let o;, c» and Cm be positive constants and let 

us consider a slightly modified (cf. the underlined term) single-dimensional primal 

formulation (5.6): Find u; € V; = H'(Qj), Ue € Ve = Hlp (Qe) such that for all 
Ve € Ve, vi € V; ° 

(c'e Vue, Vve) * Cu (Tee, T, v) C» (T;u;, TeVe) = Cn (f, Ve), 

— Cr (Teue, Tivi) + (o; Vui, Vvi) + (ui, vi) + S (Tiu Tivi) = Se (f, vi). 


(6.1) 


At At 


Here, Te, T; are the trace operators T;u; = uj|r and Tete = uelr. Letting W = Ve X V;, 
u = (ue,uj) the problem (6.1) can be stated as Au = L in W’ with the operator A 
and functional L defined as 


ee V-(aeV) + Cm TT. _ Cm Tp | T (f d | 
A = At ^e At +e and L(y) = €» > 
= TT, -V(o;V) +1 + Cu TIT; o=% =(f; "n A 


Note that the underlined lower order term in (6.1) corresponds to J in A. 
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To precondition (6.1) using the operator preconditioning framework we proceed 
to show that A : W — W’ is an isomorphism. This statement follows from the 
Lax-Milgram theorem, see e.g. (6| Ch 2.7). Let ||-|y be some norm of the space 
W (we shall see shortly that different norms can be considered leading to different 
preconditioners). Hence coercivity of A, i.e. 


There exists a, > 0 such that a. llul < (Au,u) VueW (6.3) 
and boundedness of A, i.e. 
There exists o^ > 0 such that (Au, v) < o" ||ulwllvllw. Yu,v € W (6.4) 


need to be shown. While the details of the proof are beyond the scope of the current 
text, we remark that using the space W = Hr p (Qe) x H! (0), we obtained the bound 


(Au,u) < max (ce + 2€, €x. o, 1,26; cu) llli, (6.5) 
by trace, Cauchy-Schwarz, and Young inequalities. Hence, the operator A is 
bounded. Note, however, that the constant a^ depends on Af and, in particular, 
it blows up for small time steps. Then, the lower bound is 


2 2 2,€ 2 
(Au, u) = ael Vuelly + cillVuillg + lluilo + ^ llTeue — Tiu: llo (6.6) 


. 2 
min (1,07, Ce) |lully 


and the coercivity thus holds with constant à, = min (1, @;, ce) which is indepen- 
dent of the time step. The condition number of the preconditioned system, using a 
preconditioner deduced from W, is then Z and from this theoretical consideration 
we may then expect the number of iterations to grow linearly as At decreases. 


The exact preconditioner based on W is block diagonal opearator 
_(-V(oi¥) 0 e 


5 0  -V(o,V) I 


(6.7) 


Observe that the preconditioner is efficient in the sense that its evaluation mounts to 
solving two smaller subproblems posed on Q; and Qe respectively. However, from 
our analysis we expect the performance of $ to deteriorate for small time steps. 
Indeed, this behavior is confirmed by results! in Table 6.1. Therein it can be seen 


1 All the experiments in this chapter are conducted using the manufactured problem (5.17) 
considered in the convergence study in Chapter 5. In particular, we use uniform meshes obtained 
by dividing the unit square (Qe U O;) first into n? squares each of which is then subdivided to 2 
isosceles triangles with diameter h. The finite element discretization is as described in Chapter 5. 
We set o; = 1, Ce = 2.2 and Cm = 1. 

The discrete linear systems Aj, xj, = Lp are solved iteratively using the CG or MinRes methods. 
Implementation of the methods is provided by (2). The solvers are always started from a random 
initial vector. As a convergence criterion the relative preconditioned residual norm is required to 
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that the iterations stabilize with h for At > 1074 while for smaller time steps finer 
meshes are needed to attain the bounds. 


In an attempt to find a more robust preconditioner let us observe that estimates 
(6.5), (6.6) show that u (Zu)! defines a different norm on W. Let now 
\lullw = (7t, u)!?. Then (Au,u) = lul, and (Au, v) € (Au, u)! (Av, v)? so 
that the conditions (6.3) and (6.4) hold with a* = 1, a, = 1. In another words, 
a Riesz map preconditioner with respect to the inner product (u,v)w = (lu, v) is 
independent of At. As the Riesz map in this case is in fact A! its exact evaluation 
(by LU) is not feasible. However, for the purpose of preconditioning, it suffices to 
replace the mapping by a spectrally equivalent operator. Then, provided that the 
equivalence is robust in At the approximate operator will lead to iterations bounded 
intime step and the discretization parameter. We remark that since the preconditioner 
takes the entire A into account, including the off-diagonal terms cf. block diagonal 
operator (6.7), we shall refer to it as monolithic. 


npn Diagonal (6.7) Monolithic BoomerAMG( A) 

At 2 02^ 0 pl 2 25 P NE 2? 25 2 2 
108 8 5 5 5 5 5 8 8 8 8 8 9 
106 7 7 7 7 7 7 8 8 8 8 8 9 
104 7 7 7 7 7 7 8 8 8 8 8 9 
10? 6 6 6 6 6 6 7 8 8 8 8 9 

1 Ill ua 1 TI T0 T0 8 8 8 9 9 9 
1072 26 30 33 32 33 32 8 8 9 9 10 10 
1074 39 70 90 114 139 155 | 28 26 19 14 11 10 
10-6 42 79 134 181 234 300 | 56 81 101 (104 87 62 
1078 42 78 130 179 244 327 81 140 226 320 359 378 


Table 6.1: Primal single-dimensional formulation. Number of CG iterations using 
(left) the diagonal operator (6.7) and (right) the Riesz map with respect to the A 
induced norm. Neither preconditioner is robust for At « 1. 


Since (5.6) leads to symmetric positive-definite matrices, cf. the coercivity condi- 
tion (6.3), we shall here use algebraic multigrid to construct the approximation of 
J^!. More precisely, the action of the operator is realized by single V-cycle of 
BoomerAMG (23). In Table 6.1 it can be seen that this choice leads to bounded 
CG iterations for At > 1078. However, there is a clear sensitivity of the bound to 
At. Thus BoomerAMG approximations of A“! are not At- robust. In fact, estimates 
of the condition number for the preconditioned problems with ny, = 28 are 1.2, 1.5, 
2.0, 9.2 and 53 for At = 1, 10-7, 1074, 107 and 10-8 respectively. 


be less than 107? in magnitude. Unless specified otherwise the preconditioners Bp are evaluated 
exactly by LU. Finally, condition number estimates of the preconditioned linear systems are obtained 
by using iterative Krylov-Schur solver from (10) applied to the generalized eigenvalue problem 
Apxn =An By! xy. The reported condition number is then max |4;, | / min|4;, |. 


The source code used for the experiments can be found on https: //github.com/MiroK/ 
and is arxived at (15). 
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6.2.2 Single-Dimensional Mixed Solvers 


We next consider preconditioners for the single-dimensional mixed formulation 
(5.12): Find E € S, u € Q such that 


(oJ, T) + GE (Tu J, TaT) - (U, V-T) = -(f, Tur), Vr €S, 
(q4, V- J) = 0, Vq € Q, 


where S = Hy (div, Q), Q = L? (Q) and T, is the normal trace operator TaT = T - ^i. 
Letting W = SxQ,x = (J, u) the single-dimensional mixed formulation is equivalent 
to the problem Ax = L in W’ with 


-1 At T! 


A = 0 


and L(x) = -(f,TnJ). (6.8) 
Note that in (6.8) the membrane term (T, J, TaT) is weighted by At/Cm, cf. (6.2). 
Thus, unlike in the single-dimensional primal formulation, the term does not domi- 
nate for small time steps. 


In order to apply the operator preconditioning framework to (6.8) the operator A 
needs to be shown to be an isomorphism. To this end we consider A as an abstract 
operator over W = S x Q with the form 


| (A B' h A:S—V', 69 
A = BO where BS Q., (6.9) 


and apply the Brezzi theory (7). This leads us to the potential preconditioners for the 
multi-dimensional mixed problem 


-l1;.. g al -1y7_ 2 Ata 
s= (7 I-VV ; ond s -[ I Vy tegis 0 


-1 
0 I | , (6.10) 


which are the Riesz mappings with respect to the inner products of the spaces 
W,=SxQ and W,-S$n EN xQ, (6.11) 


where N = lv € S; |v - n|lor < oo]. Note that in W2 we enforce additional regularity 
on the normal trace since T, maps S to H~!/*(L), e.g. (9] Ch. 2). Moreover the traces 
are considered in a weighted space. 


A consequence of the mapping properties of the normal trace operator is the fact that 
the term (T, J, TaT) cannot be bounded using the H (div) norm and in turn the Brezzi 
conditions are not satisfied on W;. We remark that if such a bound were possible the 
boundedness constant would depend on At, cf. (6.5). As the Brezzi conditions do not 
hold on W; we expect preconditioner 8; to perform poorly. Indeed, Table 6.2 shows 
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that the condition numbers of 8; pAn are not bounded in A. The condition numbers 
for fixed h can be seen to grow with the time step. In Table 6.3 the ill-posedness of 
the problem in W; manifests as unstable iterations or failure of MinRes to converge. 
We note that for small At the iterations appear bounded in A. 


Nh By E 
At 22 Y 27 25 P aA 27 7 27 29 26 2 
10? [4223 17540 14766 29437 4785 9569 |] 13.52 13.52 13.52 13.52 13.52 13.52 
10! 75.0 132 259 515 1028 2057]) 2.28 2.28 228 228 2.28 2.28 
1 8.86 I5 27 52 104 207 || 220 220 220 220 220 220 
107 |227 285 411 666 12 22 220 220 220 220 220 220 
107 |220 220 220 2.20 259 3.60 || 220 220 220 220 220 220 


Table 6.2: Conditioning of single-dimensional mixed formulation. (Left) Posing the 
problem in H(div) x L? violates Brezzi conditions. (Right) Preconditioner based 
on W3 in (6.11) is inf-sup stable with the inf-sup constant depending on Ar for 
At/Cm > 1. 


Posing (6.8) in W2 it can be shown that the Brezzi conditions hold with the constants 
independent of the time step. 


We shall not prove the validity of the Brezzi theory here. Instead, Table 6.2 offers nu- 
merical evidence that the discrete condition is satisfied. Therein condition numbers of 
the 82 preconditioned problems are reported and boundedness in h can be observed. 
Moreover, it can be seen that the inf-sup constant depends on Ar, in particular, the 
bound goes to 0 as At grows. However, on the subspace (Q. u) € Wo; (u, 1)o,0; = 0] 
the inf-sup condition holds independent of At. The single run-away mode appears to 
have no effect on the MinRes iterations, cf. Table 6.3, see also (20). 


np EJ EZ 

At 23 2 2 2 27 25 235 0? P 2 2 2 
108 = = - = = = 10 10 10 10 10 10 
10° - - - - = - 5 12 12 12 12 12 
104 - - B = = - 5 15 15 16 16 16 
10? 156 = = B = - 5 15 15 15 16 16 
10! 26 36 57 86 129 209 14 14 14 #14 14 14 
I 35 50 72 107 155 234 6 16 16 16 16 16 
1077 16 19 23 32 44 60 16 17 17 18 18 18 
107? 18 19 201 27 34 45 19 19 19 19 20 20 
1074 19 19 20 20 20 20 19 19 20 20 20 2 
10-6 19 19 20 20 21 21 9 19 20 20 21 2l 
1078 19 19 20 20 21 21 9 19 20 20 21 21 


Table 6.3: MinRes iterations with preconditioned single-dimensional mixed formula- 
tion using preconditioners (6.10). Failure to converge within 500 iterations is denoted 
by —. Convergence of B, seems unaffected by the single At unbounded inf-sup mode. 


Considering order optimality of the 82 based solver we recall that the iterations in 
Table 6.3 were run with the exact preconditioner. In particular, the SM N Riesz map 
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was computed by LU. In practical computations such realization is not feasible and 
scalable approximation is needed. For standard H (div) inner product, i.e. the one 
used in 84, the Riesz map can be efficiently realized by multigrid algorithms (14). 
However, in the authors' experience this approach does not work equally well with 
the S N N inner product. 


We have seen in Sections 6.2.1 and 6.2.2 that construction of order optimal monolithic 
solvers for the single-dimensional EMI formulations presents a challenge. In primal 
formulation Af robustness of the monolithic approach was problematic. For the mixed 
formulation, order optimality required a specialized solver for the Riesz mapping 
over the subspace of H (div). These two problems are addressed by solvers for multi- 
dimensional formulations. 


6.2.3 Multi-Dimensional Solvers 


In this section we consider the construction of preconditioners for the operators 


—V-(aeV) 0 -T; oI -V T; 
Ap = 0  -V(o;V) T and Am=| V 0 0 |, (6.12) 
sif; n -I T, 0 -ȘI 


which induce respectively the multi-dimensional primal weak formulation (5.8) and 
the multi-dimensional mixed weak formulation (5.14). In order to discuss well- 
posedness let 


Wp = Hip (Qe) x H' (Qi) x HPT) NYELT) (6.13) 


and 
Wm = Hys (div, Q) x L'(Q) x H^ (T) n J S (T). (6.14) 


We remark that the fractional space H^ !?, in which the membrane current density 
Im in (5.6) is sought, and H!/?, which is the space of the transmembrane potential v 
in (5.14), are dictated by the mapping properties of the trace operators Te, T; and Ty. 
For example, as T; : H'(Qj) ^ HT), Im € (H2) so that the term (7;u;, Im) is 
bounded. Note also that only the spaces involving I' are now At-dependent, cf. e.g. 
W in (6.11). 


Due to the At-weighted (penalty) terms in A, and A,, the operators cannot be 
established as isomorphisms on W, and Wm by straightforward application of the 
Brezzi theory. Instead, A, for At/C,, < 1 and Am for At/C,, 2 1 can be analyzed 
using the framework for saddle point systems with penalty, see (4| Ch. 3.4). A crucial 
assumption then is that the system without the penalty term satisfies the Brezzi 
conditions. While here the well-posedness shall be demonstrated only by numerical 
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experiments, let us point out an important difference between the operators. To this 
end let c > 0 be a constant and consider a piecewise constant potential u such that 
uj = C, Ue = Oand set J = 0. Then for any t € H(div) by divergence theorem 


(V: 7,u)o + (TaT, v) = (V-T, Wo; + (TaT, v) = C(Tnt, 1) + (TaT, v) 
and in turn, if v 2 —c we have that 
c IL-VTA(J c(T, 7,1) + (Tut, v) 0 


V. 0 Oj[u|- 0 =|0 
Tn 0 O/}\v 0 0 


Thus the operator Am without the penalty term is singular with a one dimensional 
kernel spanned by vector (0, u, —c) and the Brezzi conditions do not hold. We shall 
see implications of this property on the convergence of the iterative solvers shortly. 
We remark that A, without the penalty term is non-singular. 


Assuming that Ap : Wp — W, is an isomorphism we consider as a preconditioner 
for the single-dimensional primal formulation the operator 


~V-(aeV) 0 0 
By = 0  -V(o;V)4I 0 l (6.15) 
0 0 (CA DI? + E 


Observe that the preconditioner consists of (approximate) solvers for three decoupled 
subproblems posed on Qe, Q; and I'. Moreover, the problems on extra and intracellu- 
lar domains are standard elliptic operators for which efficient (black-box) multigrid 
techniques exist, e.g. (23). The problem on the interface is then less standard as it 
concerns a fractional Helmholtz operator. However, efficient multilevel solvers have 
been established e.g. in (5] [I). We remark that if the discrete spaces for Im (or v) 
contain only a few thousands of degrees of freedom an eigenvalue realization of the 
fractional preconditioner is feasible cf. (16). As the interfacial spaces here are small, 
cf. Table (5.1), we use further the spectral approach. 


Nh 4 5 6 7 8 39 
A: 2 03 206 0 o» 9 

tj 
ny " 3 " F n " 10 10 10 10 10 10 10 
At NE MM NE NR ME. 106 10 10 10 10 10 10 
107 2752 2752 2752 2752 2752 2752 104 I 15 15 15 15 15 
10! 29.61 29.62 29.62 29.62 29.62 29.62 102 18 18 18 17 17 17 
T 6.86 6.88 6.89 6.89 6.89 6.89 T 77 27 27 26 26 26 
1071 8.84 8.88 889 889 889 8.89 102 43 47 46 44 44 43 
10-2 | 889 8.93 894 894 894 8.94 1074 49 60 62 60 57 56 
10-6 51 60 62 61 59 59 
10-3 49 60 62 61 59 59 


Table 6.4: Fractional preconditioner (6.15) for multi-dimensional primal formulation. 
(Left) Condition numbers are bounded in A. Growth in At for At > 1 is caused by a 
single mode. (Right) MinRes iterations counts are bounded in h and time step. 
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Using (6.15) Table 6.4 reports the condition numbers of the preconditioned problems 
Bp hAp,n. It can be seen that the results are bounded in A, thus providing a numerical 
evidence for A, being an isomorphism on W, in (6.13). Moreover, it can be seen 
that the conditioning of the problem deteriorates with Aż large. However, in this 
case it is only a single mode u; = 0, ue = const, Im = 0 which causes the blowup. 
Considering the MinRes iteration counts, the presence of this mode seems to have 
almost no impact on boundedness in h or At. 


We finally come back to the multi-dimensional mixed formulation. Based on the 
space W,, in (6.14) let a preconditioner for the multi-dimensional mixed formulation 
be 
o-lI-VV-0 0 E 
Bn = 0 I 0 . (6.16) 
0 0 (CA« 1? + €x, 


We remark that the preconditioner uses a standard H (div) inner product, cf. B, in 
(6.10), which can be efficiently realized by multigrid methods, e.g. (14). Note also 
that the fractionality of the Laplacian is 1/2, cf. —1/2 of the multi-dimensional primal 
preconditioner (6.15). In addition to the previously mentioned multilevel methods, 
problems (—Au + u)?x = b for 0 < s < 1 can be efficiently solved by a number of 
approaches, see (3) and references therein. 


Using (6.16) Table 6.5 reports the condition numbers of the preconditioned multi- 
dimensional mixed formulation. The conditioning can be seen to be stable in h, 
while in agreement with the limit singularity property, there is a growth with At. 
Given that only a single mode is responsible for the lack of At-boundedness and 
recalling results of Table 6.4 or 6.3 we might expect that also here the MinRes solver 
will not be affected. However, Table 6.6 shows that this is not the case. In fact, 
as At grows the iterations become unstable in h. Figure 6.1 then shows a typical 
convergence behavior of the solver. We see that the relative preconditioned residual 
norm is quickly reduced to about 107? in approximately 30 iterations (regardless of 
the mesh resolution). Afterwards the convergence stalls. We conclude that for robust 
preconditioning the nullspace of A, must be addressed. 


nn BmAm BIAL, 

At 27 23 23 25 26 2 27 23 23 25 26 2 
107 17566 17627 17733 17810 17857 17882 || 4.40 4.42 4.51 4.58 4.63 4.67 
10! 177 178 179 180 180 180 4.39 441 450 4.57 4.62 4.66 

1 3:31 3.73 3.96 4.09 4.18 4.27 3.3] 3.73 3.96 4.09 4.18 4.27 
107! 1.08 1.16 1.30 1.56 1.96 2.48 2.61 2.61 2.61 2.61 2.64 3.14 
107? 1.03 1.03 1.03 1.03 1.03 1.03 2.62 2.62 2.62 2.62 2.62 2.62 


Table 6.5: Conditioning of the multi-dimensional mixed formulation. (Left) Problem 
considered on Wm in (6.14) yields operator Am with a singular limit as At — 
co. (Right) Formulation including an additional constraint on the transmembrane 
potential [a v dS = 0 yields h and At robust condition numbers. 
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Motivated by the observation that Am becomes singular in the limit At = co with 
a mode J = 0, ue = 0, u; = c, v = —c in the kernel, let we = We x R and 
let us define a new solution operator for the multi-dimensional mixed formulation 
JO, : W? — (W9)' and its preconditioner 89. : (W9)' > W? as 


cI-V T; 0 o-lI-VV.0 0 oy 
V 0 0 0 0 I 0 0 
Os (e 
Am=) 0 -Şa]] AmE 0 0(-A« DP + S70} " 
0 0 I1 0 0 0 0 ul 


(6.17) 
where u = min (1, Cm/At). Note that A®, includes an additional unknown, a sin- 
gle scalar, which enforces (1,v)r = 0. The new constraint thus eliminates constant 
transmembrane potential and in turn the new operator is non-singular?. Considering 
results reported in Tables 6.5, 6.6 it can be seen that the new preconditioned formu- 
lation leads to condition numbers and iteration counts, which are bounded in both 
the mesh size and the time step. Note that the extra unknown has only a small impact 
on the number of iterations compared to the Bm preconditioner. 


nn BmAm BIAS, 

At 23 23 2 oP 2 28 2j — 24-79» DD 28 
108 457 = - DA] = i - 28 32 34 36 37 37 
106 54 174 82 301 258 227 || 28. 32 34 36 37 37 
104 37 41 45 45 47 47 28 32 34 36 37 37 
10? 32 37 39 39 41 41 28 32 34 36 37 37 

T 28 30 32 34 34 35 30 33 35 35 37 38 
102 19 22 25 29 33 35 20 23 29 33 35 39 
1074 10 10 10 11 11 12 2 12 12 183 13 14 
1076 7 7 9 9 9 9 9 9 10 10 10 10 
1078 7 7 7 7 7 7 9 9 9 9 9 9 


Table 6.6: MinRes iterations for the multi-dimensional mixed formulation. No con- 
vergence in 500 iterations is indicated by —. (Left) For large At the iterations are 
unstable since operator Am becomes singular in the limit as At — co. (Right) Con- 
straining transmembrane potential the operator A?, does not have the limit singular 
property. 


? A physical motivation for the constraint can be found in considering the EMI interface equation 
Cm ge = Im + lion on E. By integrating the left-hand side we obtain Cm a, v)r. Recall that 
Im = 0; Vu; : n; on T and -V - (o; Vu;) = 0 in Q;. Then, setting Iion = 0 and integrating the 
right-hand side we have 


(1, Im + Tion)r z (1, Im)r = (1, -V. (a; Vui))o; =0. 


Thus we obtain a conservation relation Cm 4a, v)r =0. 
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Fig. 6.1: (Left) Convergence of the MinRes method for (6.12)-based multi- 
dimensional mixed formulation with Ar = 1078 and problem from Table 6.6. For 
all mesh resolutions the norm is reduced in cca. 30 iterations. Afterwards the re- 
duction stalls. (Right) Convergence of the Robin/Neumann domain decomposition 
algorithm. For small Ar the algorithm becomes sensitive to mesh size and can even 
diverge. 


6.3 Domain Decomposition Solvers 


Having seen that monolithic solvers for the EMI equations can be sensitive to spatial 
and temporal resolution we next briefly discuss robustness of the non-overlapping 
domain decomposition (DD) approach. The focus here shall only be on the single- 
dimensional primal formulation and the Robin/Neumann DD algorithm in the form 
presented in Chapter 4). In particular, our implementation shall not include any 
coarse space, cf. (21), or a preconditioner, see e.g. for interpretation of DD as 
Steklov-Poincaré operators in fractional Sobolev spaces. 


For the sake of self-containedness we review the variational formulation of the 
Robin/Neumann algorithm. Let V, = Hip (Qe), V; = H' (Qj) and v9, ul be the 
given initial transmembrane and extracellular potentials. A single iteration of the 
DD algorithm then produces a new approximation v! in three steps. (i) We find 
ui € V; such that 


(ci Vui, Vvi) + C» (Tiu; T;v;) - -Su (fv) + C» (T, ue. T; vi), Vv; € Vi. 


(ii) Using the computed intracellular potential the new extracellular potential ue € Ve 
is found such that 


(c'e Vue, Vve) = (-o;Vuj: ni Teve) Wve € Ve. 
(iii) Finally v! = T;u;- Tzu. The iterations continue by assigning v? = v! andu? = ue 


until convergence. which in the following is determined as ||v* —v*-! |lo/ ||v*-! | < e, 
where e = 1075. We remark that the tolerance was chosen such that in the refinement 
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studies, where convergence of DD was observed, the approximation error in the 
H! x H! norm decreased with the expected order as O(h) (we recall that P4 — Pj 
elements were used). 


To test robustness of the DD solver we consider the setup of the single cell experiment 
used with the monolithic approaches in Section 6.2. Here the Robin problem (i) and 
the Neumann problem (ii) shall be solved by LU to eliminate effects of inexact 
subdomain solvers. We remark that due to the conforming triangulation, see Figure 
5.1, and P, discretization of both V; and V, the discrete transmembrane potential 
is computed simply by interpolation. Figure 6.1 plots evolution of the potential 
difference DA - ue t llo/Il llo with DD iterations for several mesh resolutions and 
time step values At < 1. It can be seen that for At = 1 the algorithm converges 
in about 5 iterations irrespective of the spatial discretization. For smaller timesteps 
convergence is delayed on finer meshes and for At = 107^ the iterations diverge. 


Considering the results in Figure 6.1 we conclude that in the form presented here 
domain decomposition is not a robust algorithm for the single-dimensional primal 
formulation of the EMI equations. However, due to its simplicity (relative to e.g. 
the multi-dimensional formulations) and speed (see Section 6.4), DD might be the 
method of choice in practical cardiac modeling. It is therefore natural to ask whether 
the divergence conditions are likely to arise in real applications. To address this 
question we perform a scaling analysis based on the membrane dynamical condition 
C, 2e ~ diVuj- nj. Letting L be a characteristic length scale and V’ = LV the 
equation becomes (LCm/ r) Z = V’u; - nj and the pre-factor T = LC,,/o; can 
be seen to have the unit of seconds. Following let us insert Cm = 1uF/cm?, 
gj; = 10mS/cm, L = 100um, where the length scale is determined by the size of 
a typical cell used by the authors’. We remark that this a sensible choice given the 
setup of our experiments. As a result T = 1076s and thus At = 1 in the experiments 
reported here corresponds to a time step of 107 seconds. Note that this is the finest 
time step considered in (11). Furthermore, therein the spatial resolution is 24m 
while here with n} = 28 and L = 100um the mesh size is approximately 0.4m. In 
summary, the conditions for divergence of the DD algorithm discussed here can be 
encountered not far away from the parameter regime in (11). 


6.4 Solver Comparison 


To complete our discussion of EMI solvers we finally address the speed of the differ- 
ent algorithms. Recall that until this point results for all solvers, but the monolithic 
single-dimensional primal one, were obtained using LU in the construction. Such an 
implementation, however, is not scalable. Here we show that the proposed algorithms 
are order optimal if LU is replaced by suitable multilevel methods. 
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Referring to the legend of Figure 6.2 we shall compare 8 different methods. Single- 
dimensional primal formulation (sp) shall be solved either with the BoomerAMG(A)- 
preconditioner (mono), or diagonal preconditioner (6.7). In the latter case all the 
blocks of the diagonal preconditioner are approximated by single algebraic multigrid 
V-cycle. The single-dimensional mixed (sm) formulations shall use the 8; precondi- 
tioner (6.10) with the leading block approximated by H (div) multigrid HypreAMS 
of (14). We recall that the solver in general is not independent of the discretization, 
cf. Table 6.3, however, HypreAMS does not work well for the robust preconditioner 
Bp. Finally, the multi-dimensional primal (mp) and multi-dimensional mixed (mm) 
formulations shall use (6.15) and (6.16) with the O;, €), and Q subproblems of 
the preconditioners approximated by multigrid (BoomerAMG in case of (6.15) and 
HypreAMS for (6.16)). The fractional operators will be computed exactly. In ad- 
dition, two implementations of Robin/Neumann domain decomposition algorithm 
are considered with the subproblems solved either exactly by LU or by 4 V-cycles 
of BoomerAMG. Finally, a reference solution time is provided by timings of exact 
solution of the single-dimensional formulation (sp-LU). 


We compare the algorithms using the single cell setup of the previous experiments 
with At = 107? where this value was chosen with the intention to not favor any of 
the methods. Using tolerance e = 10? for DD and 10^? for CG/MinRes, Figure 
6.2 reports solution times? of the algorithms for 2? < np < 2!°. It can be seen that 
with the exception of LU-based solvers all the methods are indeed order optimal. 
The monolithic approach (sp-mono) and the multi-dimensional primal (mp) solver 
are then the fastest solvers while solvers for the mixed formulations (sm-5;, mm) 
are the slowest. For np = 2!? the solution times of the solvers were cca. 26s, 140s, 
640s, 570s respectively. We remark that the mixed formulations have approximately 
5 times more unknowns compared to the single-dimensional ones. Thus the cost 
per degree of freedom is comparable between the two approaches. Considering the 
scalable DD implementation (DD-AMG) the timing for n; = 2!? is cca. 160s. The 
method thus has a similar cost to multi-dimensional primal formulation. 


3 The timings include setup times of the preconditioners for monolithic methods. For DD the 
subproblems were assembled and their solvers constructed only once. 
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Fig. 6.2: Solution times (including preconditioner setup) of the solvers for the EMI 
system (5.1) in terms of number of unknowns N or mesh size A. Single cell model 
with At = 107? is used. All but multi-dimensional primal solver with LU (sp-LU) 
and LU-based domain decomposition solver (DD-LU) are order optimal. The former 
scales as N?/?. Legend is shared between the figures. 


Open Access This chapter is licensed under the terms of the Creative Com- 
mons Attribution 4.0 International License (http://creativecommons.org/licenses/ 
by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in 
any medium or format, as long as you give appropriate credit to the original aut- 
hor(s) and the source, provide a link to the Creative Commons license and indicate 
if changes were made. 

The images or other third party material in this chapter are included in the chap- 
ter's Creative Commons license, unless indicated otherwise in a credit line to the 
material. If material is not included in the chapter's Creative Commons license and 
your intended use is not permitted by statutory regulation or exceeds the permitted 
use, you will need to obtain permission directly from the copyright holder. 


References 


12. 


13. 


14. 


15. 


. Berland T, Kuchta M, Mardal KA (2019) Multigrid methods for discrete fractional Sobolev 


spaces. SIAM Journal on Scientific Computing 41(2):A948-A972 

Balay S, Abhyankar S, Adams MF, Brown J, Brune P, Buschelman K, Dalcin L, Dener A, 
Eijkhout V, Gropp WD, Karpeyev D, Kaushik D, Knepley MG, May DA, McInnes LC, Mills 
RT, Munson T, Rupp K, Sanan P, Smith BF, Zampini S, Zhang H, Zhang H (2019) PETSc Web 
page. URL 
Bonito A, Pasciak J (2015) Numerical approximation of fractional powers of elliptic operators. 
Mathematics of Computation 84(295):2083-2110 

Braess D (2007) Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics. 
Cambridge University Press, Cambridge 

Bramble J, Pasciak J, Vassilevski P (2000) Computational scales of Sobolev norms with 
application to preconditioning. Mathematics of Computation 69(230):463—480 

Brenner S, Scott R (2007) The mathematical theory of finite element methods, vol 15. Springer 
Science & Business Media 

Brezzi F (1974) On the existence, uniqueness and approximation of saddle-point problems 
arising from Lagrangian multipliers. Revue française d'automatique, informatique, recherche 
opérationnelle Analyse numérique 8(2):129-151 

Discacciati M, Quarteroni A, Valli A (2007) Robin-Robin domain decomposition methods for 
the Stokes-Darcy coupling. SIAM Journal on Numerical Analysis 45(3):1246-1268 

Girault V, Raviart PA (2012) Finite Element Methods for Navier-Stokes Equations: Theory 
and Algorithms, vol 5. Springer Berlin Heidelberg, Berlin, Heidelberg 

Hernandez V, Roman JE, Vidal V (2005) SLEPc: A scalable and flexible toolkit for the solution 
of eigenvalue problems. ACM Trans Math Software 31(3):351—362 


. Jeger KH, Tveito A (2020) Efficient numerical solution of the EMI model representing ex- 


tracellular space (E), the cell membrane (M) and the intracellular space (I) of a collection of 
cardiac cells. Preprint 

Jeger KH, Hustad KG, Cai X, Tveito A (2020) Operator splitting and finite difference schemes 
for solving the EMI model. In: Tveito A, Mardal KA, Rognes ME (eds) Modeling excitable 
tissue - The EMI framework, Simula Springer Notes in Computing, SpringerNature 

Klawonn A, Widlund OB, Dryja M (2002) Dual-primal FETI methods for three-dimensional 
elliptic problems with heterogeneous coefficients. SIAM Journal on Numerical Analysis 
40(1):159-179 

Kolev TV, Vassilevski PS (2012) Parallel auxiliary space AMG solver for H(div) problems. 
SIAM Journal on Scientific Computing 34(6):A3079-A3098 

Kuchta M, Mardal KA (2020) Software for EMI - Iterative solvers for EMI models. URL 


https://doi.org/10.5281/zenodo. 3771212 


85 


86 


16. 


17. 


18. 


19. 


20. 


21. 


22: 


23. 


24. 


25. 


Kuchta et al. 


Kuchta M, Nordaas M, Verschaeve JCG, Mortensen M, Mardal KA (2016) Preconditioners 
for saddle point systems with trace constraints coupling 2D and 1D domains. SIAM Journal 
on Scientific Computing 38(6):B962—B987 

Kuchta M, Mardal KA, Rognes ME (2020) Solving the EMI equations using finite element 
methods. In: Tveito A, Mardal KA, Rognes ME (eds) Modeling excitable tissue - The EMI 
framework, Simula Springer Notes in Computing, SpringerNature 

Mandel J, Dohrmann CR, Tezaur R (2005) An algebraic theory for primal and dual substruc- 
turing methods by constraints. Applied Nnumerical Mmathematics 54(2):167—193 

Mardal KA, Winther R (2011) Preconditioning discretizations of systems of partial differential 
equations. Numerical Linear Algebra with Applications 18(1):1—40 

Nielsen BF, Mardal KA (2013) Analysis of the minimal residual method applied to ill posed 
optimality systems. SIAM Journal on Scientific Computing 35(2):A785-A814 

Smith B, Bjorstad P, Gropp W (2004) Domain Decomposition: Parallel Multilevel Methods 
for Elliptic Partial Differential Equations. Cambridge University Press, Cambridge 

Tveito A, Jeger KH, Kuchta M, Mardal KA, Rognes ME (2017) A cell-based framework for 
numerical modeling of electrical conduction in cardiac tissue. Frontiers in Physics 5:48 

Yang UM, et al. (2002) Boomeramg: a parallel algebraic multigrid solver and preconditioner. 
Applied Numerical Mathematics 41(1):155-177 

Zampini S (2016) PCBDDC: a class of robust dual-primal methods in PETSc. SIAM Journal 
on Scientific Computing 38(5):5282-S306 

Zampini S, Tu X (2017) Multilevel balancing domain decomposition by constraints deluxe 
algorithms with adaptive coarse spaces for flow in porous media. SIAM Journal on Scientific 
Computing 39(4):A1389-A1415 


® 


Check for 
updates 


Chapter 7 


Improving Neural Simulations with the EMI 
Model 


Alessio Paolo Buccino!»?, Miroslav Kuchta?, Jakob Schreiner?, and Kent-André 
Mardal?-4 


Abstract Mathematical modeling of neurons is an essential tool to investigate neu- 
ronal activity alongside with experimental approaches. However, the conventional 
modeling framework to simulate neuronal dynamics and extracellular potentials 
makes several assumptions that might need to be revisited for some applications. 
In this chapter we apply the EMI model to investigate the ephaptic effect and the 
effect of the extracellular probes on the measured potential. Finally, we introduce 
reduced EMI models, which provide a more computationally efficient framework for 
simulating neurons with complex morphologies. 


7.1 Introduction 


In recent years, huge efforts and resources have been spent in computational mod- 
eling of neuronal activity. For example, the Blue Brain Project 
has constructed and distributed several 
hundreds of biophysically detailed cell models (multi-compartment models) from rat 
sensory cortex. A similar effort is being conducted by the Allen Institute of Brain 
Science, whose cell-type database 
includes hundreds of cell models both from mice and even from humans. As the 
experimental data used become more comprehensive and available, these models are 
expected to become elaborated and more accurate in reproducing neuronal dynam- 
ics. However, the modeling framework which is commonly used to simulate these 
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multi-compartment models makes several key assumptions that might be violated 
for certain applications. 


The most widely used approach to simulate neuronal dynamics of neurons is the 
cable equation. The solution of this equation enables one to compute transmembrane 
currents for each of the model compartment. In order to simulate extracellular 
potentials, we can use volume conduction theory and sum the individual contributions 
of the currents to the electric potential at any point in space (6). Whereas the use of 
this modeling framework has been the gold standard to simulate neuronal activity 
for decades (19}|6), there are some important assumptions that need to be discussed: 


* A neuron is represented as a discrete set of nodes. Multi-compartment models 
split the neuronal morphologies into a discrete set of segments. Therefore, neurons 
are not represented as a continuum and this might affect the accuracy of the 
simulations. However, this assumption can be alleviated by using very small 
segments that can accurately replicate the neuronal complex geometry. 


* Extracellular potentials are assumed to be constant. When solving the cable 
equation, the extracellular potentials outside the membrane are assumed to be 
constant. This assumption is harder to relax, as it prevents to include so-called 
ephaptic coupling in the simulation (9) [I). Ephaptic coupling refers to the effect 
of extracellular potentials on the neuronal dynamics. The use of the EMI model 
allows one to include these phenomena in the simulation, both to simulate the 
effect other neurons' activity or the same neuron's activity (self-ephaptic) has on 
the membrane potential. 


* The extracellular space is assumed to be homogeneous. The most common ap- 
proach to compute extracellular potentials arising from neuronal currents is to use 
volume conduction theory with the point-source or line-source approximations 
(6), which assume that the extracellular medium is homogeneous (in addition to 
linear, isotropic, and infinite). However, this assumption is clearly violated when 
using extracellular devices to record neuronal activity, which introduce a clear 
inhomogeneity in the extracellular space. Extracellular probes can be explicitly 
modeled in the extracellular space with the EMI model and they show to greatly 
contribute to the recorded signals (3). 


In applications where any of the assumptions listed above may be violated the EMI 
model (1.30)) provides a suitable modeling framework. In particular, the ge- 
ometry of the neuron (and the extracellular space) is accurately represented. Here 
we will show how the model is convenient to study both the ephaptic coupling 
of neurons (Section 7.3) and the effects of extracellular probes on the recorded 
electric potentials (Section 7.4). However, the detailed representation of the geome- 
tries makes EMI much more computationally intense than the standard modeling 
framework [3). This limits the complexity of the simulation mainly to simple 
neuronal morphologies, such as ball-and-stick neurons (3). In order to target more 
realistic morphologies Section 7.5 discusses the reduced EMI model where the 
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(three-dimensional) intracellular space is lumped to a curve that is the centerline of 
the neuron. 


7.2 EMI Simulations of Neurons using the neuronmi Python 
Package 


Before discussing the applications let us first briefly introduce a Python package 
called neuronmi!, which has been used for the simulation results reported further. 
neuronmi provides a high-level application programming interface (API) to enable 
users to easily set up and run EMI simulations of neurons. 


The workflow of the neuronmi package consists of two parts. First, a mesh needs to 
be generated. This is done with the generate_mesh function, that uses gmsh as 
backend. With this function the user can choose different kinds of neurons to place 
in the mesh and optionally place a probe device in the extracellular space. Mesh 
resolution and sizes of the bounding box can also be adjusted, as well as parameters 
of the neurons and the probe. In the following code example, we create a mesh with 
a ball-and-stick neuron (bas) and a microwire probe. By default, the center of the 
soma is at (0,0,0) um, the dendrite extends in the positive z direction and the axon in 
the negative z direction. 


import neuronmi 
mesh folder = neuronmi.generate mesh(neurons-'bas', 
probe-'microwire') 


Once a mesh is generated, the EMI simulation can be invoked with the simulate. emi 
function, which implements the finite element method for the multi-dimensional 
mixed formulation (5.14)) following the discretization proposed in (20). Through 
a set of parameters, the user can stimulate the neuron with a synaptic input, a step 
current, or a pulse current. Alternatively, the probe contacts can be used to stimulate 
the neuron extracellularly. By default, the neuron receives a synaptic input on its 
dendrite. The user can probe electric potentials u at any point in the mesh, while 
transmembrane currents i and membrane potentials v are available at facets on the 
neuron surface. The full solutions are also saved as pvd or xdmf files. The simulation 
is run as follows: 


u, i, v = neuronmi.simulate emi (mesh. folder, 
u. probe locations-points. v, 
i probe locations-points i, 
v. probe locations-points. v) 


!ihttps://github.com/MiroK/nEuronMI 
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Several parameters can be set to customize the mesh generation and the simulation. 
For further details, we refer to the package documentation (https: //neuronmi. 


readthedocs.io/en/latest/). 


7.3 Investigating the Ephaptic Effect between Neurons 


Neurons are surrounded by the electrically conductive extracellular space. Groups 
of neurons create fluctuations in the local extracellular electrical field. These fluctu- 
ations in turn influence the intracellular electrical field through the ephaptic effect. 
Ephaptic coupling cannot influence neurons at rest, however, it can affect the spike 
timings of a neuron receiving suprathreshold stimulus. 


We will illustrate how the EMI model can be used to compute the ephaptic cou- 
pling between two ball-and-stick neurons embedded in an extracellular space. The 
simulation is based on the neuronmi package detailed above. One of the neurons is 
stimulated with a synaptic input which elicits an action potential. The intracellular 
potential in the other neuron is sampled. We ran several experiments increasing the 
distance between the neurons. 


Maximum soma potential as the neurons are moved apart 


— intracellular potential 


Fig. 7.1: Two ball-and-stick neurons embedded in an extracellular space (left) and 
the increase in the intracellular potential due to ephaptic coupling (right). 


The intracellular potential is sampled in the centre of the ephaptically stimulated 
neuron (right) to measure the strength of the synaptic coupling. The deflection is 
4.7 uV when the neurons are 5 um apart and decreases to 3.2 uV when they are 40 
um apart with the soma of the stimulated neuron adjacent to the axon of the other 
neuron. 
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While in these simulations we only showed the effect of a single spike on an adja- 
cent neuron, the occurrence of synchronous activity in populations of neurons can 
cause larger degrees of ephaptic coupling. The neuronmi package enables users 
to instantiate several neurons in the mesh and to define their inputs (which can 
be synchronous), hence allowing in principle to investigate ephaptic effects at the 
population level. 


7.4 Investigating the Effect of Measuring Devices on 
Extracellular Potentials 


While the presence of recording devices is usually ignored in the computation of 
extracellular potentials, recent findings suggest that newly developed silicon-based 
devices, or Multi-Electrode Arrays (MEAs), have a strong effect on the measured 


signals (3). 


Using neuronmi, which is built on the EMI model, one can easily incorporate the 
neural devices in the mesh and investigate their effects on the recorded signals. To 
demonstrate this, we built meshes with a simple ball-and-stick neuron and different 
types of neural probes in its vicinity: 


Microwire: the first type of probe represents a microwire. For this kind of probes 
we used a cylindrical insulated model with 30 um diameter. The extracellular 
potential, after the simulations, was measured at the center of tip of the cylinder. 
The microwire probe is shown in Figure 7.2A alongside with the simplified 
neuron. 


Neuronexus (MEA): the second type of probe model represents a commercially 
available silicon MEA (A1x32-Poly3-5mm-25s-177-CM32 probe from Neu- 
ronexus Technologies), which has 32 electrodes in three columns (the central 
column has 12 recording sites and first and third columns have 10) with hexago- 
nal arrangement, a y-pitch of 18 um, and a z-pitch of 22 um. The electrode radius 
is 7.5 um. This probe has a thickness of 15 um and a maximum width of 114 um, 
and it is shown in Figure 7.2B. 


Neuropixels (MEA): the third type of probe model represents the Neuropixels 
silicon MEA (11). The original probe has more than 900 electrodes over a 
lcm shank. The probbe is 70 um wide and 20 um thick. In our mesh, shown 
in Figure 7.2C we used 24 12x12 um recording sites arranged in the chessboard 
configuration with an inter-electrode-distance of 25 um (11). 


We ran EMI simulations using the meshes with and without the probe in the extra- 
cellular space (Figure 7.2A-B-C show the meshes with the probe), and we compared 
the obtained extracellular action potentials - EAPs (Figure 7.2D-E-F, blue without 
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probe - orange with probe). The probes were placed in the extracellular space at a 
distance of 30 um from the center of the neuron. 


— no probe 
——— with probe 


da 
VY Y 
ay ye 


Fig. 7.2: Effect of different probes on the recorded potentials. (A-B-C) Meshes 
including a neuron and a microwire (A), a neuron and a Neuronexus probe, a neuron 
and a Neuropixels probe (C). (D-E-F) Extracellular action potentials computed at the 
electrodes’ location without (blue) and with the probe (orange) in the extracellular 
space. Large MEAs seem to strongly affect the recorded signals (E-F). 


Microwire probes do not affect the recorded potentials, with an EAP peak of 
—21.63 uV without the probe and of —20.53 uV with the probe (Figure 7.2D). How- 
ever, when recording with silicon MEAs, the extracellular potentials are strongly 
affected. For the Neuronexus probe (Figure 7.2D), the EAP peak without the probe 
is —30.47 uV, while with the probe it is —56.09 uV (peak ratio=1.84). For the Neu- 
ropixels probe (Figure 7.2E), the EAP peak without the probe is —32.73 uV, while 
with the probe it is —63.63 uV (peak ratio=1.94). The probe effect is probably due to 
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the insulating properties of the silicon probes, which act as electrical walls for the 
generated currents. For further details and analysis we refer to (3). 


7.5 Reduced EMI Model 


In the examples considered thus far the problem geometry was simple allowing for 
computations on a serial desktop computer. In order to apply the EMI framework to 
realistic neurons two challenges need to be addressed: representation of neurons in 
the form of a finite element mesh and efficient solvers for large linear systems due 
to the EMI equations. However, even with the order optimal algorithms discussed in 
Chapter 6) and efficient mesh generators for neuron surface geometries, see e.g. 
(17), the computational cost of the 3D-3D EMI models remains large (compared 
to the conventional approaches). As a computationally feasible alternative we shall 
next discuss the 3D-1D models. 


Topological order reduction is a modeling technique used e.g. in reservoir simula- 
tions or studies of tissue perfusion (5), which exploits geometrical properties 
of the system in order to derive its reduced model. Viewing a dendrite (branch) as 
generalized cylinder with length L and radius R we observe that R « L and that 
in a typical domain of interest the neuron’s volume is negligible compared to its 
surroundings. This property motivates a reduced representation of the neuron in 
terms of a (one dimensional) curve, the centerline, along with a function R, which 
provides radius of the crossection at each point of the line. An illustration of the 
concept can be seen in Figure 7.3. Thus, referring to the notation of Chapter 5, €2; is 
reduced to a line while Qe newly occupies the entire domain, i.e. Q = Qe. In turn, 
the reduced EMI model presents a coupling between three dimensional extracellular 
space and the one-dimensional intracellular space. We remark that the membrane is 
one-dimensional as well. 


In order to apply order reduction to the EMI model we consider the single- 
dimensional primal formulation (5.6). Note that therein, the coupling on the mem- 
brane T requires that both ue and u; are restricted from Qe and Q; respectively by 
the dedicated trace operator. In a reduced model Q; = T, T = ^ and thus u; needs 
not to be restricted. On the other hand, restriction of ue to curve A can no longer 
be realized as a trace since such an operation is not well-defined for H! functions, 
see e.g. (5). To define a value of the extracellular potential on A let us introduce an 
averaging operator 


IIu(x) = a(x) = |CaGO|[! " l O) d$, xeT,u € H!(Q). 


Here Cr(x) = (y € LT; (y - x) L n(x)} with n(x) being the unit tangent vector of A 
at x, cf. Figure 7.3. Thus, IIue is computed by sampling ue on the original (two- 
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Fig. 7.3: Reduced representation of neurons. (Left) Neuron Q; is collapsed to cen- 
terline A. Extracellular potential defined on © is evaluated at A by averaging ue 
on the cylindrical surface I'. (Right) Surface mesh and centerline representation of 
RatS 1-6-39 neuron generated by AnaMorph (17). The dendritic part satisfies R « L 
property, while the reduced model assumptions do not hold on the soma. 


dimensional) neuron surface. However, in practical computations we assume that I 
has a circular cross section so that |Cr(x)| = 2z R(x). 


Using II the reduced single-scale primal formulation of the EMI model reads: Find 
Ue € Hlp (Q), u; € H! (A) such that for all v, € Hlp (Q), vj € H'(A) 


[ove mean e f 2n ds2nR f Cu. yu; ds 
Q A A 


= ants. ds, 

A 

JE Ville ds | Roui Yvi ds+ f 2rRGeuvi ds 
A A A 


= -f 2nR fv; ds. 
A 


Here, the factors 27R arise in reducing the integration domain from I to A and 
similarly for zR? and Q;. Thus, defining the reduced specific membrane capacitance 
Cn = 2x RC, and the reduced intracellular conductivity o; = zR?o; the operator 
form of (7.1) can be seen to be (6.2) with the new restriction operators 7; = J and 
Te = II. Note also that in (7.1) the function f, which characterizes the membrane 
dynamics, is defined on I. 


(7.1) 


For the proof of well-posedness of (7.1) as well as a detailed discussion of modeling 
assumptions, which allow for model reduction from (5.6) the reader is referred to 
(4). Moreover the reduced multi-dimensional primal formulation (5.8) is analyzed 


in (13). 
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To assess the reduced model (7.1) the surface mesh of a rat neocortex neuron RatS 1-6- 
39 from the NeuroMorpho database (2) has been generated using AnaMorph (17). The 
neuron has been embedded into a (box-shaped) domain © such that |O|/|O;| ~ 100 
with the resulting geometry meshed by gmsh. The full 3D-3D single-dimensional 
primal formation has been used to compute the response to a 1 ms synaptic stimulus 
of 50 nA. Referring to Figure 7.4 the lower branch close to node number 4 has 
been stimulated. Using the centerline representation of the neuron the response has 
also been computed with (7.1). We remark that P-P; elements were used with both 
formulations and that spaces were setup on conforming meshes. In particular, with 
(7.1) the discretization of A consisted of the edges of the cells of O. However, such 
a geometrically conforming discretization is not required in the reduced model. In 
fact, the meshes of A and Q can be independent, see (4). The reduced model then 
resulted in 4587 unknowns, which is to be contrasted with 18248 unknowns due to 
(5.6). In turn, the simulation time using the reduced model is about 110 seconds 
while the full EMI model required cca. 340 seconds to complete. 


The two models are compared in Figure 7.4 which shows values of the computed 
intracellular potentials at different points along the centerline. In general, there is 
qualitative agreement between the model predictions. However, the reduced model 
can be seen to tend to underestimate both the minima and the maxima, while the 
excitation occurs faster compared to the full model. More precisely, the peak po- 
tentials due to the full model at points 2-5 were {27.64, 26.09, 22.87, 12.64} mV 
with occurrences after {2.08,2.11,2.15,2.57} ms. For the reduced model the max- 
ima {22.31,21.65, 19.39, 16.02} mV were recorded at {1.56, 1.58, 1.51, 1.86} ms. In 
addition to the intracellular potentials, the extracellular potentials were compared by 
sampling 6.52jum away from the soma center (node 1 in Figure 7.4). We remark that 
the soma radius was 5.71 jm. It can be seen that with —2.99uV < ue < 1.07uV for 
(7.1) and -1.99uV < ue < 0.67uV for (5.6) the reduced model overestimates the 
extrema. As with the extracellular potentials there is a temporal shift in the response; 
the negative peak is observed at 1.26ms, respectively 1.87ms. 


While results of the comparison in Figure 7.4 shall be viewed as preliminary we argue 
that they illustrate sufficiently the potential of reduced EMI models. In particular, 
the reduced model is able to capture qualitatively the properties of the full EMI 
simulations. However, clear differences, especially in the temporal shifts of the 
peaks, have been observed. In the future we aim to investigate if suitable scaling of 
the stimulus and/or the parameters of the membrane ODEs can reduce the prediction 
error. In addition, the modeling error of the reduced model shall be evaluated similar 
to (15). More specifically, the soma, being approximated as a sphere in the 3D-3D 
model, cannot be represented as a slender cylinder (unlike the dendrites and axons). 
Thus the model reduction assumptions are not met on the soma. While localized in 
space, the effect of this error on temporal predictions should be analyzed. In turn, 
improved reduced models, which take into account the spherical nature of the soma, 
e.g. in construction of averaging operators, might be needed. 
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Fig. 7.4: Comparison of full 3D-3D and reduced 3D-1D EMI models. RatS 1-6-39 
neuron is stimulated close to node 4 in the dendritic part (plotted in blue). Intracellular 
potentials (nodes 2-5) are measured on the centerline with node 2 being the soma (in 
red, dashed line indicates the radius) center. Extracellular potentials are compared in 
node | next to soma. Axis of each response plot is anchored at the center next to the 
measurement point. The full and reduced EMI models provide qualitatively similar 
predictions. 


7.6 Conclusions 


In this chapter, we have showcased some applications in which the EMI model 
can be a viable alternative to standard modeling frameworks in order to investigate 
aspects of neuronal activity and recorded signals. We introduced an open-source 
software package named neuronmi to easily assemble meshes including simplified 
neurons and probes in the extracellular space. Finally, we performed preliminary 
investigations into accuracy and solution cost of a reduced 3D-1D EMI model 
suitable for simulating complex neuronal morphologies. 


Open Access This chapter is licensed under the terms of the Creative Com- 
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by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in 
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